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A  Particle  Image  Velocimetry  (PIV)  system  has  been  developed  to  measure  velocity  fields 
with  order  one-micron  spatial  resolution.  The  technique  uses  200  nm  diameter  flow-tracing 
particles,  a  pulsed  Nd:YAG  laser,  an  inverted  epi-fluorescent  microscope,  and  a  cooled  interline- 
transfer  CCD  camera  to  record  high-resolution  particle-image  fields.  The  spatial  resolution  of  the 
PIV  technique  is  limited  primarily  by  the  diffraction-limited  resolution  of  the  recording  optics. 
The  accuracy  of  the  PIV  system  was  demonstrated  by  measuring  the  known  flow  field  in  a  30 
pm  x  300  pm  (nominal  dimension)  microchannel.  The  resulting  velocity  fields  have  a  spatial 
resolution,  defined  by  the  size  of  the  first  window  of  the  interrogation  spot  and  by  the  out  of 
plane  resolution  of  13.6  pm  x  0.9  pm  x  1.8  pm,  in  the  streamwise,  wall-normal,  and  out  of 
plane  directions,  respectively.  By  overlapping  the  interrogation  spots  by  50%  to  satisfy  the 
Nyquist  sampling  criterion,  a  velocity-vector  spacing  of  450  nm  in  the  wall-normal  direction  is 
achieved.  These  measurements  are  accurate  to  within  2%  full-scale  resolution,  and  are  the 
highest  spatially  resolved  PIV  measurements  published  to  date. 

In  order  to  obtain  these  measurements,  new  PIV  interrogation  algorithms  and  theory  have 
been  developed.  These  new  algorithms  include  ensemble-averaged  based  correlation  analysis 
and  central  difference  interrogation  (CDI).  The  principle  of  the  ensemble  correlations  is  that 
tremendous  gain  in  signal  to  noise  results  from  averaging  a  time-series  of  particle-image 
correlation  functions,  before  the  signal  peak  is  located.  Central  difference  interrogation  allows 
one  to  obtain  velocity  data  that  is  both  second  order  accurate  in  time  and  in  space.  The  technique 
also  allows  one  to  obtain  velocity  data  near  curved  boundaries. 

Finally,  micro  PIV  was  applied  to  obtain  measurements  of  liquid  flow  in  silicon 
micromachined  nozzles.  The  results  showed  near  separation  on  one  of  the  side  walls  in  the 
expansion  area.  Considerable  effort  was  focused  on  extending  the  technique  from  liquid  flows  in 


micronozzles  to  gaseous  flows.  This  extension  to  gaseous  flows  proved  to  be  rather  challenging. 
In  particular,  it  was  quite  difficult  to  seed  50  nm  glass  particles  in  air  flow  at  the  microscale. 
Furthermore,  imaging  these  small  particles  proved  to  be  difficult  and  non-fluorescing  imaging 
techniques,  such  as  DIC  were  explored. 

Supposed  Personnel 

Marin  Sigurdson  -  Associate  Laboratory  Engineer 
Derek  Tretheway  -  Postdoctoral  Researcher 
Steve  Wereley  -  Postdoctoral  Researcher 
Mike  Gray  -  Graduate  Student  Researcher 
Hong-Sheng  Zhang  -  Graduate  Student  Researcher 
Yang  Lu  -  Graduate  Student  Researcher 

Publications  Resulting  from  Support 

S.  T.  Wereley  and  C.  D.  Meinhart  2001.  Adaptive  second-order  accurate  particle  image 
velocimetry.  Accepted  Experiments  in  Fluids 

C.  D.  Meinhart,  S.  T.  Wereley,  and  M  .H.  B.  Gray  2000.  Volume  Illumination  for  two- 
dimensional  particle  image  velocimetry.  Measurement  Science  Technology  Volume  11,  No 
o.,  pp.  809  -  814. 

C.  D.  Meinhart,  S.  T.  Wereley,  and  J.  G.  Santiago  2000.  A  PIV  algorithm  for  estimating  time- 
averaged  velocity  fields.  Journal  of  Fluids  Engineering ,  Vol.  122,  pp.  285  -  289. 

C.  D.  Meinhart,  S.  T.  Wereley,  and  J.  G.  Santiago  1999.  PIV  Measurements  of  a  MicroChannel 
Flow.  Exp.  in  Fluids,  Vol.  27,  pp.  414-419. 

C.  D.  Meinhart,  S.  T.  Wereley,  and  J.  G.  Santiago  1999.  Micron-Resolution  Velocimetry 
Techniques.  In  Developments  in  Laser  Techniques  and  Applications  to  Fluid  Mechanics,  R. 
J.  Adrian  et  al.  (Eds.),  Springer- Verlag,  Berlin 

J.  G.  Santiago,  S.  Wereley,  C.  D.  Meinhart,  D.  J.  Beebe,  &  R.  J.  Adrian  1998.  A  PIV  system 
for  microfluidics.  Exp.  Fluids,  Vol.  25  No. 4,  pp  316-319 


Conference  Papers  Resulting  from  Support 

C.  D.  Meinhart  and  S.  T.  Wereley  2001.  Fluid  Mechanics  Issues  at  the  Microscale,  AIAA  Paper 
2001-0720,  Reno,  NV. 

S.  T.  Wereley,  L.  Gui  and  C.  D.  Meinhart  2001.  Flow  Measurements  for  the  Microfrontier, 
AIAA  Paper  2001-0243,  Reno,  NV. 

S.  T.  Wereley  and  C.  D.  Meinhart  2000.  Accuracy  Improvements  in  Particle  Image 
Velocimetry,  Algorithms  Proceedings  of  the  10th  International  Symposium  on  Applications 
of  Laser  Techniques  to  Fluid  Mechanics,  Lisbon,  Portugal. 

S.  Stone,  C.  D.  Meinhart  &  S.  T.  Wereley  2000.  A  microfluidic-based  nanoscope,  Proceedings 
of  p-TAS  2000,  Enchele,  The  Netherlands,  May  15-18 

M.  Volpert,  C.  D.  Meinhart,  I.  Mezic,  and  M.  Dahleh  1999.  An  actively  controlled  micromixer. 
ASME  -  IMECE’99  MEMS  symposium,  Nashville,  TN 

S.  T.  Wereley,  C.  D.  Meinhart,  and  M.  H.  B.  Gray  1999.  Depth  effects  in  volume  illuminated 
particle  image  velocimetry.  Proceedings  of  the  Third  International  Workshop  on  Particle 
Image  Velocimetry,  Santa  Barbara,  CA,  Sept.  16-18 

C.  D.  Meinhart,  S.  T.  Wereley,  and  J.  G.  Santiago  1999.  A  PIV  Algorithm  for  estimating  time- 
averaged  velocity  fields,  Proceedings  of  Optical  Methods  and  Image  Processing  in  Fluid 
Flow,  3rd  ASM E/JSME  Fluids  Engineering  Conference,  July  18-23,  San  Francisco,  CA 

C.  D.  Meinhart,  M.  H.  B.  Gray  &  S.  T.  Wereley  1999.  PIV  measurements  of  high-speed  flows 
in  silicon-micromachined  r  jzzles.  AIAA  Paper  No.  99-3756 

S.  T.  Wereley,  C.  D.  Meinhart,  J.  G.  Santiago  &  R.  J.  Adrian  1998.  Velocimetry  for  MEMS 
applications.  In  the  proceedings  for  the  Application  of  Microfabrication  to  Fluid  Mechanics, 
ASME  International  Mechanical  Engineering  Congress  and  Exposition,  Anaheim,  CA,  Nov. 
15-20 

J.  G.  Santiago,  S.  T.  Wereley,  C.  D.  Meinhart,  D.  Beebe,  and  R.  J.  Adrian  1998.  A  micron- 
resolution  particle  image  velocimetry  system,  8th  International  Symposium  on  Flow 
Visualization,  Sorrento,  Italy 

R.  J.  Adrian,  K.  T.  Christensen,  S.  M.  Soloff,  and  C.  D.  Meinhart  1998.  Decomposition  of 
turbulent  fields  and  visualization  of  vortices.  Proceedings  of  the  9'h  International  Symposuim 
on  Applications  of  Laser  Techniques  to  Fluid  Mechanics,  Lisbon,  Portugal,  July  13-16 


C.  D.  Meinhart,  S.  T.  Wereley,  J.  G.  Santiago,  and  R.  J.  Adrian.  1998.  Diagnostic  techniques 
for  microfluidics  research.  In  Proceedings  of  the  9,h  International  Symposium  on 
“Applications  of  Laser  Techniques  to  Fluid  Mechanics” ,  Lisbon,  Portugal,  July  13-16 
S.  T.  Werely,  J.  G.  Santiago,  R.  Chiu,  C.  D.  Meinhart  &  R.  J.  Adrian  1998  Micro-resolution 
particle  image  velocimetry.  Micro-  and  Nano-  Fabrication  Structures  and  Devices  for 
Biomedical  Environmental  Applications,  SPIE’s  BiOS’98-  International  Biomedical  Optics 
Symposium,  Jan.  24-30,  San  Jose,  CA 


1.  Micro  PIV  Development 

Recently,  there  has  been  a  growing  interest  to  develop  microscale  devices  that  can 
manipulate  and  transport  relatively  small  volumes  of  fluids.  These  devices  have  applications  in 
many  areas  of  engineering,  including  propulsion  and  power  generation  of  micro-satellites,  micro 
air  vehicles,  inkjet  printer  heads,  and  bioanalytical  instruments.  The  recent  surge  of  microfluidic 
devices  has  created  a  need  for  diagnostic  tools  with  spatial  resolutions  on  the  order  of  several 
microns.  These  diagnostic  techniques  can  be  used  as  engineering  tools  to  measure  the  flow 
performance  of  microfluidic  devices,  and  as  general  scientific  instruments  to  probe  the  physics  of 
transport  processes  at  the  microscale. 

Particle  image  velocimetry  (PIV)  is  a  well-established  technique  for  measuring  velocity 
fields  in  macroscopic  fluid  systems  (Adrian,  1991).  In  double-pulsed  PIV,  positions  of  flow¬ 
tracing  particles  are  recorded  at  two  known  times  by  either  illuminating  the  particles  using  a 
pulsed  light  source,  or  by  illuminating  the  particles  using  a  continuous  light  source  and  gating  the 
light  near  the  camera  using  a  mechanical  or  electronic  shutter.  The  displacement  of  the  particle 
images  is  then  estimated  statistically  by  correlating  the  particle  image  pairs  (Meinhart  et  al, 
1993). 

The  first  successful  micro  PIV  experiment  was  conducted  by  Santiago  et  al.,  (1998).  In  that 
experiment,  we  used  an  epi-fluorescent  microscope  with  a  continuous  Hg-arc  lamp,  and  a 
Princeton  Instruments’  intensified  CCD  camera  to  record  the  flow  around  a  nominally  30  pm 
cylinder  in  a  Hele-Shaw  flow  cell.  A  bulk  velocity  of  50  pm  s'1  was  measured  with  a  spatial 
resolution  of  6.9  pm  x  6.9  pm  x  1.5  pm,  based  upon  the  size  of  the  first  correlation  window  and 
the  depth-of-field  of  the  microscope.  The  flow-tracing  particles  were  chosen  to  be  300  nm 
diameter  polystyrene  particles.  These  particles  were  large  enough  to  scatter  sufficient  light  for 
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recording  and  to  reduce  the  effects  of  Brownian  motion.  The  Hg-arc  lamp  continuously 
illuminated  the  particles,  and  the  light  signal  was  gated  electronically  by  the  CCD  intensifier 
plate.  In  this  experiment,  the  CCD  array  was  exposed  for  bt  =  2  ms,  and  the  time  delay  between 
exposures  was  At  =  68.5  ms. 

The  experimental  setup  of  Santiago  et  al.  (1998)  is  well  suited  for  situations  where  high 
spatial  resolution  and  low-light  levels  are  required,  such  as  investigating  flows  around  living 
microorganisms.  Unfortunately,  the  minimum  exposure  time  of  the  CCD  camera  necessary  to 
record  particle  images  is  on  the  order  of  several  milliseconds,  and  the  time  delay  between  image 
exposures  is  on  the  order  of  tens  of  milliseconds.  This  large  time  interval  between  image 
exposures  limits  the  PIV  system  to  relatively  low  velocities. 

In  this  paper  we  describe  an  alternative  approach  where  a  5  ns  pulsed  Nd:YAG  laser  is  used 
to  illuminate  sub-micron  fluorescent  particles  with  ~1  mJ  of  light  energy.  The  fluorescent  images 
are  recorded  using  a  cooled  interline  transfer  1300  x  1030  x  12  bit  CCD  camera  that  is  capable 
of  taking  back-to-back  images  within  a  time  interval  as  short  as  500  ns.  The  CCD  camera’s  pixel 
spacing  is  6.8  Jim. 

Flow-Tracing  Particles 

Developing  PIV  to  achieve  microscale  spatial  resolution  requires  that  the  particles  are  chosen 
small  enough  to  follow  the  flow  faithfully  without  (1)  disrupting  the  flow  field,  (2)  clogging  the 
microdevice,  and  (3)  producing  unnecessarily  large  images.  At  the  same  time,  the  particles  must 
be  chosen  large  enough  so  that  they  scatter  sufficient  light  to  be  recorded  and  sufficiently 
dampen  out  Brownian  motion. 

Making  accurate  PIV  measurements  with  spatial  resolutions  on  the  order  of  several  microns 
requires  that  the  diameter  of  flow-tracing  particles,  dp,  be  on  the  order  of  100  nm  to  300  nm 


diameter.  If  visible  light  with  a  wavelength  of  X  =  532  nm  is  used  to  illuminate  the  particles,  the 
particle  diameter  will  be  smaller  than  the  wavelength  of  light,  i.e.  dp<\.  In  this  regime,  it  is 
difficult  to  image  particles  using  normal  elastic  scattering  techniques.  However,  inelastic 
scattering  techniques,  such  as  epi-fluorescence,  can  be  used  to  image  sub-micron  particles. 
Fluorcscently  labeled  polystyrene  particles  with  diameters  of  100  -  300  nm  and  a  specific  gravity 
of  pp=  1.055  are  suitable  for  many  microfluidic  liquid-flow  applications.  The  particles  can  be 
doped  with  a  fluorescent  dye  that  is  tuned  to  absorb  light  from  a  frequency-doubled  Nd:YAG 
laser  (A  =  532  nm),  with  a  peak  excitation  wavelength  of  Xabs  =  540  nm  and  a  peak  emission 
wavelength  of  Xemit  =  560  nm.  When  sub-micron  particles  are  used  to  trace  slow  flows,  one  must 
consider  errors  due  to  particle  diffusion  resulting  from  Brownian  motion.  A  first  order  estimate 
of  this  error  relative  to  the  displacement  in  the  x-direction  is  given  by  Santiago  et  al.  (1998) 


u 


(1) 


where  s2  is  the  random  mean  square  particle  displacement  associated  with  Brownian  motion, 
and  D  is  the  Brownian  diffusion  coefficient,  u  is  the  characteristic  velocity,  and  At  is  the  time 
interval  between  pulses. 

Since  the  errors  due  to  Brownian  motion  are  unbiased,  they  can  be  substantially  reduced  by 
averaging  over  several  particle  images  in  a  single  interrogation  spot  or  by  ensemble  averaging 
over  several  realizations.  In  practice,  errors  due  to  Brownian  motion  place  a  lower  limit  on  the 
size  of  particle  that  can  be  used  to  achieve  the  desired  velocity  measurement  accuracy.  This  size 
is  dependent  on  the  characteristic  velocity  of  the  experiment. 

The  required  density  of  flow-tracing  particles  in  micro  PIV  experiments  is  commonly  on  the 
same  order  as  that  for  macro  PIV  experiments.  Consider  a  macro  PIV  experiment  where  one  uses 


a  50  |im  diameter  particle  to  seed  a  water  flow,  and  an  interrogation  spot  size  of  1  mm  x  1  mm  x 
1  mm.  Assuming  an  average  of  10  particles  per  interrogation  volume  are  required  for  adequate 
measurement  signal  (using  standard  autocorrelation  techniques),  the  ratio  of  particle  volume  to 
interrogation  spot  volume  is  ~  0.07%.  In  the  current  experiment,  there  was  an  average  of  2.5 
particles  per  a  13.6  pm  x  0.9  pm  x  1.8  pm  interrogation  volume.  Given  particle  diameters  of  200 
nm,  the  ratio  of  particle  volume  to  interrogation  spot  volume  is  0.05%.  In  order  to  obtain  high 
quality  signal  with  only  2.5  particles  per  interrogation  spot,  we  employed  a  specialized 
ensemble-averaged  correlation  technique.  The  algorithm  is  used  to  average  over  20  instantaneous 
correlation  functions,  before  searching  for  the  signal  peak.  This  algorithm  increases  the  number 
of  effective  particles  per  measurement  volume,  and  provides  adequate  signal  for  reliable 
measurements.  A  complete  description  of  the  ensemble-averaged  correlation  technique  is 
discussed  by  Meinhart  et  al.  (1999). 

Imaging  and  Recording 

All  optical  diagnostic  techniques  are  ultimately  limited  by  the  diffraction  of  light.  For  micro 
PIV,  the  spatial  resolution  is  limited  by  the  effective  diameter  of  particle  images  when  projected 
back  into  the  flow  field.  For  magnifications  much  larger  than  unity,  the  diameter  of  the 
diffraction-limited  point  spread  function,  in  the  image  plane,  is  given  by 

d  =  2.44  A/ ,  (2) 

1  2  NA 

where  M  is  the  total  magnification  of  the  microscope  and  NA  is  the  numerical  aperture  of  the 
lens  (Bom  &  Wolf,  1997).  Using  an  oil-immersion  Nikon  CFI  Plan  Apochromat  60x  objective 
lens  with  a  numerical  aperture  of  NA  =  1.4,  and  assuming  the  wavelength  of  the  recording  light 
is  X  ~  560  nm,  the  diameter  of  the  point  spread  function  in  our  experiment  is  ds  ~  29.3  pm. 
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The  actual  image  recorded  on  the  CCD  camera  is  the  convolution  of  the  diffraction-limited 
image  with  the  geometric  image  (Adrian,  1991).  Approximating  both  the  geometric  and 
diffraction-limited  images  as  Gaussian  functions,  the  resulting  convolution  is  a  Gaussian 
function  with  an  effective  particle  diameter  de,  where 

de=[M'dl+dlf .  (3) 

For  a  magnification  of  M  =  60,  a  particle  diameter  of  dp  =  200  nm,  and  numerical  aperture  of 
NA  =  1.4,  the  effective  particle  image  diameter  projected  onto  the  CCD  camera  is  de  =  31.7  pm. 
The  effective  particle  diameter  when  projected  back  into  the  flow  is  528  nm.  According  to 
Prasad  et  al.  (1993)  if  a  particle  image  diameter  is  resolved  by  3  -  4  pixels,  the  location  of  a 
particle-image  correlation  peak  can  be  determined  to  within  l/10th  the  particle-image  diameter. 
This  yields  a  measurement  uncertainty  of  5x  =  dJIOM.  For  the  parameters  consider  here,  the 

uncertainty  reduces  to  Sx  ~  528 nm  / 10  =  52.8  nm. 

The  fact  that  one  can  measure  particle  displacement  to  within  53  nm  is  somewhat  surprising. 
In  most  microscopic  applications,  one  is  primarily  interested  in  determining  the  shape  of  small 
objects.  Obviously,  the  smallest  resolvable  shape  is  on  the  order  of  the  resolution  of  the 
microscope.  In  micro  PIV,  one  knows  a  priori  the  particle  shape  and  is  interested  only  in 
determining  particle  position.  By  over-sampling  the  image  (i.e.  resolving  the  image  with  3-4 
pixels  across  the  image  diameter),  one  can  determine  particle  position  to  within  an  order  of 
magnitude  better  resolution  than  the  diffraction-limited  resolution  of  the  microscope. 

In  most  PIV  applications,  the  thickness  of  a  laser  light  sheet,  A z,  determines  the  out-of-plane 
measurement  domain.  This  thickness  is  usually  chosen  to  be  smaller  than  the  depth-of-field  of 
the  recording  system,  5 z.  Consequently,  all  particles  illuminated  by  the  light  sheet  produce  in¬ 
focus  images,  reducing  background  noise  in  the  image  field  (Adrian,  1991).  In  micro  PIV,  one  is 
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often  interested  in  obtaining  velocity  measurements  with  out-of-plane  resolutions  on  the  order  of 
1-10  /tm.  Forming  a  light  sheet  in  a  microchannel  that  is  1  -  10  /tm  thick  is  extremely  difficult. 
Aligning  the  light  sheet  so  that  it  overlaps  with  the  in-focus  object  plane  of  the  image-recording 
system  is  nearly  impossible.  In  principle,  a  light  sheet  could  be  formed  in  a  microfluidic  device 
by  micromachining  optical  wave  guides  and  light  sheet  forming  optics  directly  into  the  device 
(Sobek  et  al.,  1994).  Because  of  the  added  difficulty  of  integrating  micromachined  optics  into 
our  system,  we  have  chosen  an  alternative  method  that  uses  broad  field  illumination. 

Figure  1  shows  schematic  of  the  micro  PIV  system.  The  illumination  beam  is  produced  by 
two  Nd:YAG  lasers  manufactured  by  New  Wave,  Inc.  The  beam  is  adjusted  by  a  set  of  lenses  and 
delivered  into  an  inverted  epi-fluorescent  microscope,  where  an  optical  filter  assembly  directs 
the  beam  to  the  objective  lens.  The  objective  lens  relays  the  light  onto  the  microfluidic  device, 
where  it  illuminates  the  entire  flow  volume.  Fluorescent,  particles  in  the  flow  field  absorb  the 
green  illumination  light,  X  =  532  nm,  and  emit  a  distribution  of  red  light  (approximately  X  ~  560 
nm).  The  emitted  light  is  imaged  through  an  oil-immersion  lens  and  passed  to  the  fluorescent 
filter  cube,  where  green  light  from  background  reflections  is  filtered  out  and  the  ~ed  fluorescence 
from  the  sub-micron  particles  is  recorded  onto  a  Princeton  Instruments  cooled  interline-transfer 
CCD  camera.  The  camera  is  cooled  to  T  =  -15°C ,  which  lowers  the  readout  noise  so  the  weak 
fluorescent  signal  can  be  measured.  The  duration  of  the  laser  pulse  is  5  ns.  Within  500  ns  after 
exposure,  the  image  field  is  transferred  to  storage  pixels  on  the  CCD  camera,  so  that  a  second 
image  field  can  be  recorded  by  the  CCD  camera.  After  a  specified  time  delay,  At,  a  second  laser 
pulse  is  used  to  record  a  second  set  of  particle  images  onto  the  CCD  camera.  Both  of  these 
images  are  then  read  out  of  the  CCD  camera  and  downloaded  to  a  PC  computer  for  processing. 
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The  depth  of  field  of  a  microscope  system  is  somewhat  arbitrary,  because  there  are  several 
definitions  of  when  an  image  is  considered  to  be  unfocused.  Inoue  and  Spring  (1997)  estimate 


the  total  depth  of  field  as  the  sum  of  the  depth  of  field  due  to  diffraction  and  geometric  affects 

(4) 
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where  n  is  the  index  of  refraction  of  the  immersion  medium  between  the  microfluidic  device  and 
the  objective  lens  (n  =  1.515  for  an  oil  immersion  lens),  X  is  the  wavelength  of  light  in  a  vacuum 
(k  ~  560  nm  for  our  fluorescent  dye),  NA  is  the  numerical  aperture  of  the  objective  lens,  M  is  the 
total  magnification  of  the  system,  and  e  is  the  smallest  resolvable  distance  of  the  image  detector 
(i.e.  the  pixel  spacing  for  the  case  of  a  CCD  camera). 

The  first  term  in  equation  (4)  is  due  to  diffraction.  The  length  represented  by  this  term  is 
chosen  by  convention  to  be  one  fourth  of  the  total  distance  between  the  first  intensity  minima 
along  the  optical  axis  of  the  3-D  point  spread  function  (Inoue  &  Spring,  1997).  The  NA-2 
dependence  of  this  term  suggests  that  the  depth  of  field  for  a  high  NA  imaging  system  decreases 
more  rapidly  than  the  lateral  resolution  of  the  image,  as  NA  increases.  For  the  present 
experiment,  equation  (4)  yields  8z~  0.6/rm . 

Equation  (4)  should  not  be  directly  applied  to  determine  out  of  plane  resolution  for  micro 
PIV  systems.  A  better  criterion  for  out  of  plane  resolution  is  to  determine  the  distance  along  the 
optical  axis  when  a  particle  becomes  sufficiently  unfocused  so  that  it  only  contributes  a  fraction, 
say  1/1 0th,  to  the  correlation  function,  compared  to  a  similar  particle  that  is  located  at  the  object 
plane.  In  our  experiment,  we  estimated  the  out  of  plane  resolution  by  focusing  the  objective  lens 
on  a  set  of  particles  fixed  to  a  microscope  slide,  and  recording  a  series  of  images  with  the 
objective  lens  placed  at  different  axial  positions.  The  out  of  plane  resolution  was  then  estimated 
by  determining  the  distance  the  object  plane  moved  to  produce  sufficiently  out  of  focus  particle 
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images,  which  could  not  significantly  contribute  to  particle-image  correlation.  We  estimated  our 
out  of  plane  resolution  to  be  Az  «  1.8  ±  0.5 /J.  m . 

The  process  to  determine  out  of  plane  resolution  should  be  conducted  in  a  timely  fashion  to 
minimize  drift  in  the  microscope  (Crenshaw,  1998).  Since  the  human  eye  can  accommodate  from 
about  25  cm  to  infinity,  out  of  plane  resolution  and  depth  of  field  cannot  be  determined 
accurately  from  direct  visual  observation  (Inoue  &  Spring,  1997). 

Lastly,  out  of  focus  particle  images  produce  background  noise  in  the  image  field,  which  can 
decrease  the  signal  of  the  in-focus  particles.  In  order  to  obtain  reliable  measurements,  the 
background  noise  must  be  minimized  by  using  high  dynamic  range  CCD  cameras,  applying 
special  image  processing  algorithms  during  interrogation,  and  limiting  the  number  density  of 
particles  in  the  flow  field. 

MicroChannel  Measurements 

A  30  pm  x  300  pm  x  25  mm  glass  rectangular  microchannel,  fabricated  by  Wilmad 
Industries,  was  mounted  flush  to  a  170  pm  thick  glass  coverslip  and  a  microscope  slide  (see 
Figum  2).  By  carefully  rotating  the  glass  coverslip  and  the  CCD  camera,  the  channel  was 
oriented  to  the  optical  plane  of  the  microscope  within  0.2°,  in  all  three  angles.  The  orientation 
was  confirmed  optically  by  focusing  the  CCD  camera  on  the  microchannel  walls.  The 
microchannel  was  horizontally  positioned  using  a  high-precision  x-y  stage,  and  verified  optically 
to  within  ~  400  nm  using  epi-fluorescent  imaging  and  image  enhancement. 

The  glass  microchannel  was  imaged  through  an  inverted  epi-fluorescent  microscope  and  a 
Nikon  Plan  Apochromat  oil-immersion  objective  lens  with  a  magnification  M  =  60  and  a 
numerical  aperture  NA  =  1.4.  The  object  plane  was  placed  at  approximately  7.5  ±  lpm  from  the 
bottom  of  the  30  pm  thick  microchannel.  The  Plan  Apochromat  lens  was  chosen  for  the 
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experiment,  because  it  is  a  high  quality  microscope  objective  designed  with  low  curvature  of 
field,  low  distortion,  and  corrected  for  spherical  and  chromatic  aberrations  (Inoue  &  Spring, 
1997). 

Since  deionized  water  was  used  as  the  working  fluid,  the  effective  numerical  aperture  of  the 
objective  lens  was  limited  to  NA  *  1.23 .  A  filtered  continuous  white  light  source  was  used  to 
align  the  test  section  with  the  CCD  camera  and  to  test  for  proper  particle  concentration.  During 
the  experiment,  the  continuous  light  source  was  replaced  by  the  pulsed  Nd:YAG  laser.  A 
Harvard  Apparatus  syringe  pump  was  used  to  produce  a  200/jI  hr~'  flow  through  the 
microchannel. 

The  particle-image  fields  were  analyzed  using  a  PIV  interrogation  program  developed 
specifically  for  microfluidic  applications.  The  program  uses  an  ensemble-averaging  correlation 
technique  to  estimate  velocity  vectors  at  a  single  measurement  point  by  (1)  cross  correlating 
particle-image  fields  from  twenty  instantaneous  realizations,  (2)  ensemble  averaging  the  cross 
correlation  functions,  and  (3)  determining  the  peak  of  the  ensemble-averaged  correlation 
function.  This  process  is  repeated  for  each  velocity  vector  in  the  measurement  domain.  The 
signal-to-noise  ratio  is  significantly  increased  by  ensemble  averaging  the  correlation  function 
before  peak  detection,  as  opposed  to  either  ensemble  averaging  the  velocity  vectors  after  peak 
detection,  or  ensemble  averaging  the  particle-image  field  before  correlation.  The  ensemble¬ 
averaging  correlation  technique  is  limited  to  steady  or  periodic  flows.  For  the  current 
experiment,  twenty  realizations  were  chosen  because  that  was  more  than  a  sufficient  number  of 
realizations  to  give  excellent  signal,  even  with  a  first  interrogation  window  of  only  120  x  8 
pixels. 
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The  signal-to-noise  ratio  resulting  from  the  ensemble-average  correlation  technique  was  high 
enough  that  there  were  no  erroneous  velocity  measurements.  Consequently,  no  vector  validation 
was  performed  on  the  data  after  interrogation.  The  velocity  field  was  smoothed  using  a  3  x  3 
Gaussian  kernel  with  a  standard  deviation  of  1  grid  spacing  in  both  directions. 

Figure  3  shows  an  ensemble-averaged  velocity- vector  field  of  the  microchannel.  The  images 
were  analyzed  using  a  low  spatial  resolution  away  from  the  wall,  where  the  velocity  gradient  is 
low,  and  using  a  high  spatial  resolution  near  the  wall,  where  the  wall-normal  velocity  gradient  is 
largest.  The  interrogation  spots  were  chosen  to  be  longer  in  the  streamwise  direction  than  in  the 
wall-normal  direction.  This  allowed  for  a  sufficient  number  of  particle  images  to  be  captured  in 
an  interrogation  spot,  while  providing  the  maximum  possible  spatial  resolution  in  the  wall- 
normal  direction.  The  spatial  resolution,  defined  by  the  size  of  the  first  interrogation  window  was 
120  x  40  pixels  in  the  region  far  from  the  wall,  and  120  x  8  pixels  near  the  wall.  This 
corresponds  to  a  spatial  resolution  of  13.6  pm  x  4.4  pm  and  13.6  pm  x  0.9  pm,  respectively.  The 
interrogation  spots  were  overlapped  by  50%  to  satisfy  the  Nyquist  sampling  criterion. 
Consequently,  the  velocity-vector  spacing  in  the  wall-normal  direction  was  450  ran  near  the 
wall.  The  streamwise  velocity  profile  was  estimated  by  line-averaging  the  measured  velocity 
data  in  the  streamwise  direction.  Figure  4  compares  the  streamwise  velocity  profile  estimated 
from  the  PIV  measurements  (shown  as  symbols)  to  the  analytical  solution  for  laminar  Newtonian 
flow  in  a  rectangular  channel  (shown  as  a  solid  line).  The  agreement  is  within  2%  full-scale 
resolution.  The  bulk  flow  rate  of  the  analytical  curve  was  determined  by  matching  the  free- 
stream  velocity  data  away  from  the  wall.  The  wall  position  of  the  analytical  curve  was 
determined  by  extrapolating  the  velocity  profile  to  zero  near  the  wall. 
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Since  the  microchannel  flow  was  fully  developed,  the  wall-normal  component  of  the  velocity 
vectors  should  be  close  to  zero.  The  average  angle  of  inclination  of  the  velocity  field  was  found 
to  be  small,  0.0046  radians,  suggesting  that  test  section  was  slightly  rotated  relative  to  the  CCD 
array.  This  was  corrected  mathematically  by  rotating  the  coordinate  system  of  the  velocity  field 
by  0.0046  radians.  The  precise  location  of  the  wall  was  more  accurately  determined  by  applying 
the  no-slip  boundary  condition  and  extrapolating  the  velocity  profile  to  zero  at  16  different 
streamwise  positions  (see  Figure  3a). 

Most  PIV  experiments  have  difficulty  measuring  velocity  vectors  very  close  to  the  wall.  In 
many  situations,  hydrodynamic  interactions  between  the  particles  and  the  wall  prevent  the 
particles  from  traveling  close  to  the  wall,  or  background  reflections  from  the  wall  overshadow 
particle  images.  By  using  200  nm  diameter  particles  and  epi-fluorescence  to  remove  background 
reflections,  we  have  been  able  to  make  accurate  velocity  measurements  to  within  about  450  nm 
of  the  wall,  see  Figure  4.  The  position  of  the  wall  can  be  determined  to  within  400  nm  by  direct 
observation  of  the  image. 
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2.  PIV  Algorithm  Theory 


Ensemble- Average  Correlation  Analysis 

At  the  microscale,  it  is  difficult  to  form  a  light  sheet  that  is  only  a  few  microns  thick, 
and  even  more  difficult  to  align  such  a  light  sheet  with  the  objective  plane  of  the 
recording  optics.  Therefore,  in  micro  PIV  experiments  the  entire  test  section  volume  is 
illuminated  with  a  cone  of  light  emanating  from  the  recording  lens,  instead  of  a 
commonly  used  light  sheet.  This  limits  the  number  density  of  particles  that  can  be  used  to 
trace  the  flow.  If  the  particle  density  is  too  high,  then  background  noise  from  out-of-focus 
particles  will  dominate  the  image  and  reduce  the  visibility  of  in-focus  particles.  If  the 
particle  density  is  too  low,  then  standard  cross  correlation  techniques  will  fail  to  provide 
an  adequate  signal,  causing  the  measurements  to  be  noisy  and  sometimes  unreliable. 

Since  many  of  the  low  Reynolds  number  flows  in  microfluidic  devices  are  laminar 
and  either  steady  or  periodic,  it  is  not  necessary  to  determine  instantaneous  velocity 
information.  For  many  of  these  flows,  it  is  sufficient  to  measure  only  the  time-averaged 
or  phase-averaged  velocity  field.  In  principle,  these  flow  fields  can  be  measured  using 
high-resolution  pointwise  techniques  such  as  the  dual  beam  laser  Doppler  anemometer 
(LDA)  system  developed  by  Tieu,  Machenzie  &  Li  (1995).  This  LDA  system  uses  low  f- 
number  optics  to  focus  the  probe  volume  down  to  approximately  5  pm  x  10  pm.  The 
micron-resolution  LDA  technique  has  advantages  over  the  micro-PIV  technique  in  that 
one  does  not  have  to  deal  with  out  of  focus  particle  images,  and  that  one  can  obtain  large 
numbers  of  velocity  measurements  in  real  time.  Micro  PIV  has  advantages  over  micro 
LDA  in  that  it  does  not  require  alignment  of  laser  beams  inside  a  microfluidic  device,  and 
it  does  require  sweeping  the  probe  volume  throughout  the  measurement  domain.  In 


17 


addition,  micro  PIV  takes  advantage  of  the  high  signal  to  noise  offered  by  fluorescence 
imaging. 

We  present  a  PIV  algorithm  that  directly  estimates  time-averaged  or  phase-averaged 
velocity  fields.  The  algorithm  can  provide  reliable  measurements  in  low  signal-to-noise 
conditions,  where  standard  cross  correlation  techniques  fail.  The  PIV  algorithm  is 
demonstrated  by  measuring  flow  in  a  30  pm  x  300  pm  microchannel  flow. 

Decomposition  of  the  Correlation  Function 

The  auto  correlation  function  of  a  single-frame  double-exposure  particle  image  field, 
I(X ),  is  defined  as 

^(s)  =  JJ/(X)/(A'  +  s)d2X  ,  (5) 

where  X  is  the  spatial  coordinate  in  the  image  plane  and  s  is  the  spatial  coordinate  in  the 
correlation  plane. 

R(s )  can  be  decomposed  into  the  following  components 

R(s)  =  Rc  (s )  +  Rp  (s)  +  Rt  (s)  +  Ro_  (s)  +  Rp+  (s) ,  (6) 

where  Rc  is  the  convolution  of  mean  image  intensity,  Rp  is  the  pedestal  component 
resulting  from  each  particle  image  correlating  with  itself,  and  RF  is  the  fluctuating 
component  resulting  from  the  correlation  between  fluctuating  image  intensity  with  mean 
image  intensity  (Adrian,  1988).  The  positive  and  negative  displacement  components  are 
R  and  Rd_  ,  respectively. 

The  signal  used  to  measure  displacement  is  contained  only  in  the  R^  component. 
Consequently,  contributions  from  all  other  components  of  the  correlation  function  may 
bias  the  measurements,  add  random  noise  to  the  measurements,  or  even  cause  erroneous 
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measurements.  Therefore,  it  is  desirable  to  reduce  or  eliminate  all  the  other  components. 
If  both  exposures  of  a  double-pulse  particle  image  field  are  recorded  separately  as  a 
double-frame  image  field,  then  the  particle  images  can  be  interrogated  using  cross 
correlation  (Keane  &  Adrian,  1992).  In  the  case  of  double-frame  cross  correlation,  the  Rp 
and  R  components  do  not  ^xist.  Further  improvements  in  signal  quality  can  be  obtained 
by  image  processing.  By  subtracting  the  mean  image  intensity  at  each  interrogation  spot 
before  correlation,  the  Rc  and  Rf(s)  components  are  set  identically  to  zero,  leaving  only 
the  positive  displacement  component,  ,  in  the  correlation  function. 

While  the  R  .  component  is  commonly  considered  the  signal  component,  it 

contains  both  noise  and  signal  information.  In  order  to  make  accurate  measurements,  the 
signal-to-noise  ratio  in  the  R  component  must  be  sufficiently  high.  This  paper  presents  a 
method  in  which  the  signal-to-noise  ratio  can  be  maximized  when  one  is  interested  in 
measuring  average  displacement  fields. 

Methods  of  Estimating  Average  Velocity  Fields 

Average  velocity  fielus  can  be  obtained  by  first  measuring  instantaneous  velocities 
and  then  averaging  them  in  either  space  or  time.  Here  we  present  an  alternative  method 
for  measuring  average  velocity  fields  directly  from  the  correlation  function. 

Estimation  of  velocity-vector  fields  using  PIV  involves  three  primary  steps: 

1 .  Particle  Image  Acquisition 

2.  Particle  Image  Correlation 

3.  Correlation  Peak  Detection 
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In  order  to  obtain  an  average  velocity  measurement,  one  must  apply  an  averaging 
operation.  The  averaging  operator  is  a  linear  operator,  and  can  be  applied  after  step  (1), 
step  (2),  or  step  (3),  to  produce  a  non-biased  estimate  of  average  velocity.  The  particle- 
image  correlation  and  peak  detection  operations  are  both  nonlinear,  and  the  order  in 
which  the  averaging  operator  is  applied  can  dramatically  change  the  quality  of  the 
resulting  signal. 

In  this  section  we  examine  the  three  different  methods  for  calculating  average 
velocities  based  upon  the  order  in  which  the  averaging  operator  is  applied.  Here  we 
assume  that  two  single  exposure  images,  Image  A  and  Image  B,  are  separated  by  a  known 
time  delay,  At,  and  represent  the  image  acquisition  of  a  single  realization.  Furthermore,  a 
sequence  of  two  images  is  acquired  at  N  different  realizations  of  a  statistically  stationary 
process.  We  wish  to  average  over  the  N  realizations  to  estimate  the  average  velocity. 

Average  Velocity  Method 

In  this  method,  the  estimate  of  the  average  velocity  is  determined  by,  (1)  correlating 
Image  A  and  Image  B,  (2)  detecting  the  peak  from  the  instantaneous  correlation 
functions,  and  (3)  averaging  over  the  instantaneous  velocity  measurements.  Figure  5a 
graphically  depicts  this  process.  The  primary  advantage  of  this  method  is  that  one  can 
obtain  instantaneous  velocity  measurements,  which  may  be  of  physical  importance. 
However,  if  one  is  primarily  interested  in  the  average  velocity  field,  this  method  may  not 
be  optimal. 

The  nonlinear  operation  of  peak  detection  is  susceptible  to  producing  erroneous 
measurements  when  the  signal  to  noise  is  low  in  the  instantaneous  correlation.  In 
practice,  the  vector  fields  must  be  validated  by  identifying  and  removing  erroneous 


velocity  measurements  (Meinhart  et  al.,  1994;  Westerweel,  1993).  Without  vector 
validation,  all  the  instantaneous  measurements  must  be  reliable  in  order  to  obtain  a 
reliable  average  velocity  measurement.  For  situations  where  the  particle  image  density  is 
low,  there  may  not  be  adequate  signal  to  obtain  valid  measurements  from  the 
instantaneous  correlation  functions.  In  these  situations,  the  alternative  methods  described 
below  may  be  employed. 

Average  Image  Method 

In  this  method,  the  averaging  operation  can  be  performed  directly  on  image  fields  A 
and  B,  and  then  correlated  to  obtain  R-B .  This  process  is  depicted  graphically  in  Figure 
5b.  In  the  situation  where  particle  image  number  density  is  low,  averaging  the  image 
fields  together  can  increase  the  average  number  of  particles  per  interrogation  spot. 
However,  excessive  averaging  of  the  image  fields  can  create  too  many  particle  images  in 
the  interrogation  spot  and  reduce  the  visibility  of  individual  particle  images. 

The  2-D  spatial  correlation  function  of  the  average  of  A  and  B  image  fields  is  defined 
as 

«JS(s)  =  }Jl(X)5(X  +  s)dJX,  (7) 

_  J  N 

where  A  =  —  ^  Ai  is  the  average  operator,  averaging  over  N  realizations  of  the  flow 
field.  Expanding  the  integrand  of  (7)  yields 
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^(‘)=lf 

Following  the  discussion  of  previous  section,  if  A  and  B  represent  paired  single¬ 
exposure  image  fields  separated  by  a  known  time  interval,  At,  which  are  cross  correlated, 
and  if  the  mean  image  intensity  is  subtracted  before  correlation,  then  all  NxN  terms  in 
(8)  are  contained  in  the  Rr>,  component  of  R. 

Only  the  N  diagonal  terms  in  (8)  contribute  to  signal.  The  N{N-l)  off  diagonal  terms 
represent  random  particle  correlation  across  different  realizations,  which  produce  noise  in 
the  correlation  function  that  can  reduce  the  accuracy  of  the  measurement  and  cause 
erroneous  measurements.  Since  the  off  diagonal  terms  are  random  and  uniformly 
distributed  throughout  the  correlation  plane,  they  may  or  may  not  contaminate  the  signal 
peak.  The  number  of  noise  terms  increases  as  n 2 ,  while  the  number  of  signal  terms 
increases  as  N. 

In  demanding  situations  when  particle-number  density  is  low,  or  when  small 
interrogation  regions  are  chosen  to  produce  high  spatial  resolution,  it  may  be  impossible 
to  obtain  an  adequate  signal  from  instantaneous  cross  correlation  analysis.  Clearly, 
averaging  the  particle-image  fields  over  several  realizations  and  then  correlating  can 
increase  the  signal-to-noise  ratio.  However,  if  the  number  of  realizations  used  in  the 
average  is  beyond  an  optimal  number,  the  signal-to-noise  ratio  can  be  reduced  by  further 
increases  in  N.  Loss  of  signal  occurs  because  the  number  of  non-diagonal  (noise)  terms 


A\BX  +  A{B2  +  AXB2  +•••+  A^Bf;  + 
A2Bi  +  A1B1+A1Bl  +•••+  A2BN  + 

AiBl  +  A2B2  +  A2B3  +•••+  A2Bn  + 

•  •  • 

•  •  • 

•  •  • 

anb{ + anb2  +  anb3  +•  •  •+ anbn 


d'X  . 


(8) 
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increases  faster  than  the  diagonal  (signal)  terms.  The  optimal  number  of  realizations 
depends  upon  many  factors  including  interrogation  spot  size,  particle-image  number 
density,  and  particle-image  quality. 

A  verage  Correlation  Method 

A  third  method  for  estimating  average  displacement  fields  is  to  calculate 
instantaneous  correlation  functions,  average  the  correlation  function,  and  then  determine 
the  location  of  the  signal  peak  location.  The  Average  Correlation  Method  is  shown 
graphically  in  Figure  5c.  The  average  of  instantaneous  correlation  functions  over  N 
realizations  can  be  written  as 

R„(s)~tfA(X)B(X  +  S)d-X 
=  j]zl(X)5(A:  +  s)d2X 

Figure  6  shows  several  instantaneous  correlation  functions,  and  an  average 
correlation  function  calculated  by  averaging  20  instantaneous  correlation  functions.  The 
instantaneous  correlation  functions  are  characteristic  of  micro  PIV  data  and  were 
obtained  from  an  actual  micro  PIV  experiment.  The  instantaneous  correlation  functions 
contain  significant  amounts  of  noise  that  can  lead  to  inaccurate  or  unreliable 
measurements.  This  noise  is  substantially  less  in  the  average  correlation  function. 

Since  the  operations  of  averaging  and  the  integration  commute,  the  average  operator 
can  be  taken  inside  of  the  integral  so  that 

Rab(s)  =\][A^B^+A2B1+AiB^+»»»  +  ANBN'\&1  X  .  (10) 

From  (10),  the  average  correlation  function  contains  only  the  diagonal  (i.e.  signal) 
terms  of  (8).  Therefore,  this  type  of  average  produces  a  much  higher  signal-to-noise  ratio 
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than  the  alternative  methods  that  have  been  discussed  previously.  The  number  of  signal 
terms  increase  linearly  with  the  number  of  realizations,  N .  Unlike  the  average  image 
method,  there  is  no  penalty  for  including  large  numbers  of  realizations  in  the  average. 
Since  the  average  operator  is  applied  to  the  correlation  function  before  peak  detection,  the 
probability  of  erroneous  measurements  is  greatly  reduced. 

Application  to  MicroChannel  Flow 

The  three  different  averaging  algorithms  were  applied  to  a  series  of  particle  image 
fields  taken  of  steady  Stokes’  water  flow  through  a  30  p.m  x  300  (im  glass  microchannel. 
Details  of  the  same  experiment  but  for  a  different  flow  rate  are  discussed  by  Meinhart  et 
al.  (in  press).  The  out  of  plane  measurement  resolution  is  related  to  the  depth  of  field  of 
the  microscope  optics,  and  is  estimated  to  be  approximately  1.8  (im.  Based  on  the 
discussion  of  Santiago  et  al.  (1998),  the  error  due  to  Brownian  motion  for  a  single 
particle  is  estimated  to  be  less  than  3%  full-scale  error.  The  interrogation  spot  size  was  16 
x  64  pixels  and  24  x  72  pixels,  for  window  1  and  window  2,  respectively.  The  relative 
offset  between  the  windows,  was  adjusted  adaptively  to  approximate  the  particle  image 
displacement  at  each  measurement  location.  The  interrogation  windows  were  chosen  to 
be  small  and  have  a  high  aspect  ratio  to  achieve  high  spatial  resolution-  in  the  wall-normal 
direction. 

The  signal-to-noise  ratio  from  standard  correlation  techniques  was  relatively  low, 
because  there  was  an  average  of  only  2.5  particle  images  located  in  each  16  x  64  pixel 
interrogation  window.  The  signal-to-noise  ratio  could  be  improved  by  increasing  the  size 
of  the  correlation  windows,  but  with  reduced  spatial  resolution.  Alternatively,  the  particle 
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concentration  in  the  flow  field  could  be  increased,  but  that  would  increase  the  number  of 
out  of  focus  particles,  and  produce  more  background  image  noise. 

Figure  7a  shows  an  instantaneous  velocity-vector  field,  calculated  using  a  standard 
cross  correlation  technique  without  any  type  of  vector  validation  or  smoothing  applied  to 
the  field.  The  velocity  measurements  contained  in  the  row  closest  to  the  wall  are  not 
considered  to  be  accurate,  because  part  of  their  interrogation  region  lies  outside  the  flow 
domain.  The  velocity  measurements  are  noisy,  and  approximately  20%  appear  to  be 
erroneous.  An  average  over  20  instantaneous  velocity  fields  is  shown  in  Figure  7b. 
Averaging  reduces  small  random  noise  in  the  instantaneous  fields.  However,  the  large 
errors  associated  with  the  erroneous  measurements  are  not  averaged  out  quickly.  The 
application  of  vector  validation  algorithms  can  usually  eliminate  errors  due  to  erroneous 
measurements  producing  accurate  velocity  measurements,  however  a  priori  knowledge 
of  the  flow  field  is  required. 

The  same  series  of  double-frame  particle  image  fields  were  correlated,  the  correlation 
functions  were  averaged,  and  then  the  locations  of  the  signal  peaks  in  the  correlation 
plane  were  determined,  following  the  average  '■'Or  relation  method  depicted  in  Figure  5c. 
The  resulting  velocity-vector  field  is  shown  in  Figure  7c.  No  vector  validation  or 
smoothing  has  been  applied  to  this  field.  The  field  contains  only  0.5  -  1  %  erroneous 
measurements  and  has  less  noise  than  the  vector  fields  shown  in  Figures  7a  and  7b. 

The  relative  performance  of  the  three  averaging  algorithms  was  quantitatively 
compared  by  varying  the  number  of  realizations  used  in  the  average,  from  one  to  twenty. 
The  fraction  of  valid  measurements  for  each  average  was  determined  by  identifying  the 
number  of  velocity  measurements  in  which  the  streamwise  component  of  velocity 
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deviated  by  more  than  10%  from  the  known  solution  at  each  point.  For  this  comparison, 
the  known  solution  was  the  velocity  vector  field  estimated  by  applying  the  average 
correlation  technique  to  20  realizations,  and  then  smoothing  the  flow  field. 

Figure  8  shows  the  fraction  of  valid  measurements  for  each  of  the  three  algorithms  as 
a  function  of  the  number  of  realizations  used  in  the  average.  The  average  correlation 
method  outperforms  the  other  two  methods,  and  produces  less  than  0.5  -  1%  erroneous 
measurements  after  averaging  eight  realizations.  The  average  image  method  produces 
about  95%  reliable  velocity  measurements,  and  reaches  a  maximum  at  four  averages. 
Further  increases  in  the  number  of  realizations  used  to  average  the  images  decreases  the 
signal  to  noise  of  the  average  particle-image  field,  and  produces  noise  in  the  correlation 
plane  due  to  random  correlation  between  non-paired  particle  images.  The  average 
velocity  method  reaches  a  maximum  of  88%  reliable  measurements  using  two  velocity 
averages.  Further  increases  in  the  number  of  averages  decreases  the  number  of  reliable 
measurements,  due  to  an  increase  in  the  probability  of  an  encountering  an  erroneous 
measurement. 
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Central  Difference  Interrogation 


Since  the  early  days  of  particle  image  velocimetry  (PIV),  researchers  have  realized  that  PIV 
algorithms  do  not  generally  make  maximal  use  of  the  information  contained  in  the  images  they 
interrogate.  In  particular,  Spedding  and  Rignot  (1993)  observed  that  the  most  appropriate  spatial 
location  for  the  displacement  vector  calculated  by  measurement  techniques  which  track  the 
motion  of  tracer  particles,  such  as  PIV,  is  displaced  from  the  centroid  of  the  first  interrogation 
region  by  a  distance  equal  to  half  the  displacement  vector,  and  not  at  the  centroid  of  the  first 
interrogation  region  (as  is  commonly  implemented).  Attributing  the  displacement  measurement 
to  its  proper  location  is  equivalent  to  a  central  difference  approximation  (order  At2)  to  the 
velocity  field  whereas  the  conventional  implementation  of  PIV  is  equivalent  to  a  forward 
difference  approximation  (order  At).  Wereley,  et  al,  (1998)  and  Lourenco  and  Krothapalli, 
(1998)  both  implemented  second  order  accurate  PIV  techniques  without  proof  of  their 
effectiveness.  Unfortunately,  moving  the  measured  displacement  vector  from  the  centroid  of  the 
first  interrogation  region  destroys  the  uniform  spacing  of  the  measured  displacement  data  that 
most  researchers  prefer.  Because  other  velocimetric  techniques,  such  as  particle  tracking 
velocimetry  and  scalar  image  velocimetry,  also  produce  non-uniformly  spaced  data,  several 
researchers  in  fields  related  to  PIV  have  addressed  how  to  map  non-uniformly  spaced  data  onto  a 
uniform  grid  and  what  are  the  effects  of  doing  so  (Spedding  and  Rignot,  1993,  Agui  and 
Jimenez,  1987,  Cohn  and  Koochesfahani,  1999).  Regardless  of  the  accuracy,  the  remapping 
process  will  always  introduce  some  error  into  the  measurement. 

In  this  paper  we  will  explore  another  approach,  namely  adaptively  shifting  the 
locations  of  the  first  and  second  interrogation  regions  such  that  when  the  calculated 
displacement  vectors  are  located  at  their  most  proper  locations  (displaced  from  the 
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centroids  of  the  first  interrogation  regions  by  half  the  displacement  distance),  they  will 
occur  on  a  uniformly  spaced  grid,  preserving  the  accuracy  of  the  original  calculations  and 
eliminating  the  need  for  remapping.  We  will  present  an  adaptive,  iterative,  second  order 
accurate  PIV  algorithm  which  significantly  improves  spatial  resolution  and  overall 
accuracy  while  reducing  velocity  bias  errors  in  comparison  to  conventional  PIV 
algorithms.  This  algorithm  was  initially  developed  by  Wereley,  et  al.  (1998)  and 
presented  without  demonstrating  its  second  order  accuracy.  Here  the  analysis  is 
expanded  and  a  theoretical,  computational  and  experimental  analysis  of  the  algorithm 
will  be  given.  We  will  consider  two  component  PIV  in  which  the  images  are  recorded  as 
pairs  of  singly-pulsed  images  on  consecutive  frames  of  some  digital  medium  (CCD 
camera  typically).  Most  of  the  ideas  considered  here  are  extensible  to  other  forms  of  two 
component  and  three  component  PIV. 

Spatial  Resolution 

The  spatial  resolution  of  any  physical  measuring  system  can  be  generally  defined  as  the  volume 
jf  space  required  to  make  the  measurement.  Similarly,  the  temporal  resolution  of  a  measurement 
system  can  be  defined  as  the  total  amount  of  time  over  which  the  measurement  is  averaged. 

Using  these  definitions,  the  spatial  resolution  of  cross-correlation  PIV  is  determined  by  the  size 
of  the  interrogation  windows  used  to  interrogate  the  images  and  the  spatial  shift  between  the 
interrogation  windows.  The  temporal  resolution  is  the  time  between  each  image  exposure. 
Generally,  the  spatial  resolution  K  is  a  three-dimensional  function  of  position  that  can  be 
expressed  as 

{kx  ,  Ky ,  Kz  }=  {min(Wxl ,  Wx2 )+  |AX|,  min(wyl ,  Wy2 )+  |AY|,  5zc ) 
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where  WxJ  and  Wy)  represent  the  size  of  the  first  interrogation  window  and  Wx2  and  Wy2 
represent  the  size  of  the  second  interrogation  window.  The  spatial  shift  between  first  and  second 
interrogation  windows  in  the  x  direction  is  represented  by  AX  and  in  the  y  direction  by  AY .  The 
quantity  8zc  represents  the  depth  of  the  measurement  domain  in  the  out  of  plane  direction.  These 
parameters  can  be  constants  and  express  the  overall  spatial  resolution  for  the  entire  image  or  they 
can  vary  from  measurement  point  to  measurement  point  and  produce  a  spatial  resolution  which 

4 

varies  across  the  region  of  interest. 

The  typical  constraint  that  the  difference  between  the  particle  displacement  and  the  spatial  offset 
between  first  and  second  interrogation  regions  not  exceed  1/3  of  the  interrogation  region  size  can 
be  eliminated  by  always  using  a  spatial  shift  equal  to  the  local  flow  velocity.  The  window  size 
W  can  be  chosen  by  signal  strength  or  spatial  gradient  considerations  alone.  Keane  and  Adrian 
(1993)  were  the  first  to  suggest  such  an  approach  and  performed  the  theoretical  calculations  to 
show  that  this  method  would  be  superior  to  other  methods.  Urushihara,  Meinhart,  and  Adrian 
(1993)  first  implemented  a  simple,  flexible  spatial  shifting  technique.  The  second  order  accurate 
adaptive  technique  presented  here  is  similar  in  nature  to  the  first  order  accurate  technique  of 
Cowen  and  Monismith  (1997),  who  implement  a  multi-pass  adaptive  spatial  shift  as  part  of  their 
super  resolution  particle  tracking  routine.  At  every  measurement  point  they  started  with  an 
initial  estimate  of  the  velocity  at  that  point  which  was  used  as  the  spatial  shift.  The  resulting 
velocity  measurement  was  then  used  in  the  next  iteration  as  the  spatial  shift.  Measurement 
points  at  which  this  procedure  did  not  converge  to  a  constant  spatial  shift  within  three  iterations 
were  flagged  as  invalid. 
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Accuracy 


The  issue  of  accuracy  is  related  to  that  of  spatial  resolution.  Larger  interrogation  windows  and 
larger  spatial  shifts  will  clearly  produce  more  coarse  estimations  of  the  local  or  instantaneous 
velocity  field.  However,  simply  minimizing  the  ci  -e  of  the  interrogation  windows  or  the  distance 
the  particles  move  between  the  two  images  is  not  the  solution  because  a  certain  minimum 
number  of  particle  images  are  required  to  produce  reliable  measurements  in  any  given 
experimental  situation.  In  addition,  those  particles  must  move  a  sufficient  distance  to  reduce  the 
relative  error  in  determining  the  location  of  the  correlation  peak  to  acceptable  levels.  One 
solution  is  to  use  a  higher-order  approximation  velocity  calculation.  To  take  advantage  of  the 
existing  specialized  equipment  for  acquiring  image  pairs — double  pulse  lasers  and  interline 
transfer  cameras — it  is  necessary  to  limit  the  data  acquisition  to  pairs  of  images.  Most  PIV 
algorithms  use  a  simple  forward  differencing  interrogation  (FDI)  scheme  in  which  the  velocity  at 
time  t  is  calculated  using  particle  images  recorded  at  times  t  and  t+  At,  given  by 


V(t)  = 


dX 

dr 


X(t  + At)-X(f)  A t  d^X 
A t  ■  2  dz2 


+  . 


(12) 


The  accuracy  of  the  velocity  calculation  can  be  greatly  improved  by  using  a  relatively  simple 
central  differencing  interrogation  (CDI)  scheme  in  which  the  velocity  at  time  t  is  calculated  using 
particle  images  at  times  t-At/2  and  t+  At  /2,  given  by 


"(')-f 


X(t  + At/2)-X(t- At/2)  (At)  d^X 
At  24  dz 3 


+  . 


(13) 


The  velocity  approximation  is  accurate  to  order  At  for  the  forward  difference  formulation  while 
the  central  difference  formulation  is  accurate  to  order  At2  using  the  same  two  images.  This 
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improved  accuracy  can  correspond  to  better  accuracy  for  a  given  time  delay  or  allow  increasing 
the  time  delay  between  images  while  maintaining  a  desired  accuracy.  The  relative  contribution 
of  other  errors,  such  as  Brownian  motion  or  the  uncertainty  in  finding  the  correlation  peak 
location,  to  the  overall  error  can  then  be  reduced. 

Velocity  Bias 

One  area  where  measurement  errors  can  be  significantly  reduced  at  little  additional 
computational  expense  is  velocity  bias.  PIV  measurements  made  with  equally-sized 
interrogation  windows  in  the  first  and  second  images  are  generally  biased  towards  the  spatial 
shift  of  the  interrogation  windows.  Westerweel  (1997)  demonstrates  that  this  error  is  on  the 
order  of  0. 1  pixels  for  common  PIV  parameters — 32x32  pixel  interrogation  windows,  particle 
diameters  in  the  range  of  2  to  4  pixels,  10  particle  images  per  interrogation  window,  and 
displacements  ranging  from  0  to  10  pixels.  Westerweel  also  demonstrates  that  this  error  decays 
linearly  to  zero  as  the  particle  displacement  decreases  from  0.5  pixels  to  zero  pixels.  While  this 
error  may  seem  relatively  insignificant,  it  can  seriously  degrade  the  calculation  of  derivative 
quantities  such  as  vorticity  and  shear  stress.  Even  though  the  velocity  bi<  error  becomes  very 
small  for  displacements  smaller  than  0.5  pixels,  Huang,  et  al.  (1997)  demonstrate  that  simply 
spatial  shifting  identically-sized  interrogation  windows  is  insufficient  to  fully  eliminate  velocity 
bias.  The  velocity  bias  can  be  eliminated  by  several  techniques.  Westerweel  (1993) 
demonstrates  that  dividing  the  image  correlation  function  by  the  correlation  of  the  windowing 
functions,  commonly  denoted  Fi,  will  completely  eliminate  the  velocity  bias.  Recently,  Huang, 
et  al.  (1997)  demonstrated  an  efficient  algorithm  that  eliminates  velocity  bias  errors  by 
renormalizing  the  values  of  the  correlation  function  in  the  vicinity  of  the  peak  location  prior  to 
calculating  the  subpixel  peak  location. 
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Keane  and  Adrian  (1993)  suggest  that  bias  errors  can  be  nearly  eliminated  by  making  the  second 
interrogation  window  sufficiently  larger  than  the  first  so  that  all  particle  images  in  all  first 
interrogation  windows  will  likely  be  contained  in  the  second.  They  also  propose  that  spatially- 
shifting  the  second  interrogation  window  by  the  integer  part  of  the  displacement  will 
substantially  reduce  this  error.  In  a  further  analysis  of  Keane  and  Adrian’s  theory,  Westerweel 
(1997)  points  out  that  using  interrogation  windows  of  differing  sizes  results  in  a  constant-valued 
portion  at  the  center  of  the  window  correlation  function  Fi,  corresponding  to  small 
displacements,  which  eliminates  the  velocity  bias.  Since  our  adaptive  procedure  uses  spatial 
shifts  that  very  accurately  reflect  the  velocity  of  the  particles  in  the  interrogation  windows,  the 
approach  suggested  by  Keane  and  Adrian  can  be  implemented  with  minimal  computational 
expense  by  making  one  of  the  interrogation  windows  as  little  as  2  pixels  larger  than  the  other  in 
each  direction.  Because  the  application  of  our  technique  is  measuring  microscale  flows  in  which 
Brownian  motion  of  the  seed  particles  can  be  a  serious  factor,  we  typically  make  the  second 
interrogation  window  larger  than  the  first  to  account  for  possible  diffusion  of  the  tracer  particles. 
The  size  of  the  second  interrogation  window  is  increased  by  the  number  of  pixels  that  particle 
diffusion  and  spatial  gradients  are  expected  to  affect  the  flow  plus  a  few  additional  pixels, 
typically  two,  in  each  direction  to  produce  the  constant-valued  region  at  the  center  of  the  window 
correlation  function. 

Ensemble-Averaged  Correlation  Function 

For  both  the  experimental  and  simulated  flows  presented  here,  an  ensemble  of  image  pairs  is 
acquired  or  generated  which  can  be  interrogated  at  higher  spatial  resolution  than  single  pairs  of 
images.  Consequently  the  effect  of  the  spatial  averaging  inherent  in  PIV  measurements  is 
reduced,  allowing  the  accuracy  of  the  central  differencing  algorithm  to  be  fully  explored.  The 
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standard  method  for  determining  ensemble-averaged  flow  fields  is  to  obtain  a  series  of 
instantaneous  velocity  fields  and  then  ensemble  average  them.  When  the  signal-to-noise  ratio  of 
the  images  is  low,  which  can  happen  for  a  variety  of  reasons,  such  as  faint  seed  particle  images, 
high  background  noise,  and  Brownian  motion  of  seed  particles,  this  procedure  can  produce  many 
spurious  vectors,  called  ‘outliers’,  due  to  the  peak  search  algorithm  mistakenly  identifying  a 
noise  peak  as  the  signal  peak.  These  outliers  can  be  located  anywhere  in  the  correlation  plane 
and  are  completely  uncorrelated  with  the  signal  peak  location.  Consequently  they  can  cause 
significant  noise.  Often  these  vectors  can  be  eliminated  from  the  ensemble  by  using  statistical 
arguments.  Even  so,  it  is  preferable  for  these  outliers  not  to  be  created  in  the  first  place. 

Recently,  Meinhart,  Wereley,  and  Santiago  (2000)  demonstrated  a  novel  ensemble  averaging 
technique  for  improving  the  signal-to-noise  ratio  of  steady-state  velocity  measurements.  Instead 
of  ensemble  averaging  after  the  calculation  of  the  series  of  velocity  fields,  at  every  measurement 
point  they  ensemble  averaged  the  correlation  function  before  the  peak  search  algorithm  was 
performed.  By  ensemble  averaging  in  this  way,  the  signal  peaks  add  constructively  while  the 
randomly  distributed  noise  peaks  do  not,  making  the  signal  peak  more  distinguishable  from  the 
noise  peaks.  This  procedure  computes  an  ensemble  averaged  velocity  field  without  first 
calculating  the  instantaneous  velocity  field.  One  of  the  primary  advantages  of  ensemble 
averaging  in  this  manner  is  that  the  effective  particle  image  concentration  is  increased,  which 
allows  for  reducing  the  size  of  the  interrogation  windows.  In  the  experiments  and  simulations 
described  in  Section  4,  26  image  pairs  were  acquired  or  generated  which  allowed  the  spatial 
resolution  to  be  increased  to  the  point  at  which  20  statistically-independent  measurements  can  be 
made  across  the  relatively  narrow  gap  of  a  cylindrical  Couette  flow  system. 
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Central  Difference  Interrogation  Fundamentals 


A  central  difference  interrogation  (CDI)  scheme  has  many  advantages  over  the  conventional 
forward  difference  interrogation  (FDI)  scheme.  The  most  obvious  advantage  already  mentioned 
is  that  the  technique  is  second  order  accurate,  providing  increased  accuracy  over  the  forward 
difference  scheme  for  the  same  time  delay  or,  alternately,  the  same  accuracy  as  the  forward 
difference  scheme  for  a  significantly  longer  time  delay.  In  many  situations  it  is  advantageous  to 
use  a  larger  time  delay.  For  example,  in  microscale  flows  the  diffusion  of  submicron  seed 
particles  may  be  a  significant  source  of  error  for  microscale  velocity  measurements  (Santiago,  et 
al.,  1998).  The  error  associated  with  Brownian  motion  eB  was  shown  to  be 
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where  u  is  the  local  velocity,  D  is  the  diffusion  coefficient  of  the  particles,  and  At  is  the  time 
delay.  Equation  14  demonstrates  that  when  all  other  factors  remain  constant  the  relative  error 
due  to  Brownian  motion  can  be  reduced  by  increasing  the  time  delay. 


(14) 


Further  advantages  of  the  CDI  algorithm  are  that  it  allows  velocity  field  to  be  calculated  closer  to 
flow  boundaries  as  well  as  in  a  symmetric  manner  with  respect  to  flow  boundaries.  These 
advantages  can  be  demonstrated  by  examining  the  following  stagnation  type  flow.  Figures  9a 
and  9b  show  velocity  fields  of  a  flow  around  a  human  red  blood  cell  calculated  by  an  adaptive 
spatial  shift  FDI  scheme  and  also  by  the  adaptive  spatial  shift  CDI  scheme.  These  results  were 
taken  from  a  set  of  experiments  conducted  by  Juan  Santiago  at  the  Beckman  Institute  at  the 
University  of  and  are  described  further  by  Wereley,  et  al.  (1998).  The  blood  cell  is  trapped 
between  a  coverslip  and  a  slide  glass  and  surrounded  by  deionized,  particle  seeded  water  (300 
nm  diameter  polystyrene  latex,  10"4  concentration  by  volume).  The  flow  is  driven  from  lower 
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left  to  upper  right  around  the  blood  cell  using  the  surface  tension  of  a  water  droplet  placed  at  the 
edge  of  the  coverslip  to  establish  a  flow  past  the  blood  cell.  The  images  are  interrogated  with 
28x28  pixel  interrogation  windows  that  are  overlapped  by  75%.  While  this  overlap  oversamples 
the  images,  it  functions  as  a  very  accurate  (yet  computationally  inefficient)  interpolator  of  the 
velocity  field.  For  both  algorithms,  measurement  regions  that  resulted  in  more  than  20%  of  the 
combined  area  of  the  first  and  second  interrogation  windows  inside  of  the  blood  cell  are 
eliminated  and  replaced  with  an  ‘x’  symbol  because  they  will  tend  to  produce  velocity 
measurements  with  serious  errors.  The  CDI  scheme  is  able  to  accurately  measure  velocities 
closer  to  the  surface  of  the  cell  than  is  the  FDI  scheme.  The  CDI  scheme  has  a  total  of  55  invalid 
measurement  points,  equivalent  to  59.4  pm2  of  image  area  that  cannot  be  interrogated,  while  the 
forward  difference  has  57  invalid  measurement  points  or  61.6  pm2.  Although  this  difference  of 
two  measurement  points  may  not  seem  significant,  it  amounts  to  the  area  that  the  FDI  algorithm 
cannot  interrogate  being  3.7%  larger  than  the  area  the  CDI  algorithm  cannot  interrogate. 
Furthermore,  the  distribution  of  the  invalid  points  is  significant.  By  carefully  comparing  Figure 
9a  with  9b,  it  is  apparent  that  the  FDI  has  three  more  invalid  measurement  points  upstream  of  the 
blood  cell  than  does  the  CDI  scheme  while  the  CDI  scheme  has  one  more  invalid  point 
downstream  of  the  blood  cell.  This  difference  in  distribution  of  invalid  points  translates  into  the 
centroid  of  the  invalid  area  being  nearly  twice  as  far  from  the  center  of  the  blood  cell  in  the  FDI 
case  (0.66  pm  or  7.85%  of  the  cell  diameter)  versus  the  CDI  case  (0.34  pm  or  4.05%  of  the  cell 
diameter).  This  difference  means  that  the  FDI  measurements  are  less  symmetrically  distributed 
around  the  blood  cell  than  the  CDI  measurements  are.  In  fact,  they  are  biased  toward  the  time 
the  first  image  was  recorded.  Computing  the  average  distance  between  the  invalid  measurement 
points  and  the  surface  of  the  blood  cell  indicates  how  closely  to  the  blood  cell  surface  each 
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algorithm  will  allow  the  images  to  be  accurately  interrogated.  On  average  the  invalid 
measurement  points  bordering  the  red  blood  cell  generated  using  the  FDI  scheme  are  12%  farther 
from  the  cell  than  the  invalid  measurement  points  generated  with  the  CDI  scheme. 

Consequently,  the  adaptive  CDI  algorithm  is  more  symmetric  than  the  adaptive  FDI  algorithm 
and  also  allows  measuring  the  velocity  field  nearer  the  cell  surface. 

Central  Difference  Interrogation  Algorithm 

The  forward  difference  scheme  (conventional  PIV)  can  be  easily  implemented  using  a  constant 
spatial  shift  in  a  single  iteration,  while  the  central  difference  scheme  requires  a  more  complicated 
algorithm.  In  general,  the  goal  of  PIV  is  to  measure  a  velocity  field  that  contains  spatial 
gradients.  These  spatial  gradients  translate  into  particles  moving  at  different  velocities  as  a 
function  of  spatial  position.  In  order  to  define  the  measurement  time  t  as  halfway  between  t-At/2 
and  t+At/2,  the  spatial  shift  between  interrogation  windows  at  each  measurement  point  in  the  first 
and  second  image  must  be  equal  to  the  distance  the  particles  move  between  images,  just  as  in 
conventional  forward  difference  algorithms.  In  practice,  a  uniformly-spaced  grid  of 
measurement  points  is  defined  at  the  time  t  and  the  interrogation  windows  in  the  first  image  are 
shifted  by  an  amount  -  V(x,  y)  •  At  12  while  the  interrogation  windows  in  the  second  image  are 
shifted  by  an  amount  +  V(x,  y)  •  At  /2,  netting  a  spatial  shift  of  V(x, y)  •  At  centered  between  the 
particle  locations  in  the  first  and  second  images.  Clearly,  employing  this  variable  spatial  shift 
requires  an  advanced  knowledge  of  the  velocity  field.  An  iterative  adaptive  approach  is  required 
in  which  increasingly  accurate  estimates  of  the  velocity  field  are  used  to  compute  the  spatial 
shift,  in  turn  yielding  increasingly  accurate  velocity  field  estimates. 
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This  procedure  is  outlined  in  flowchart  form  in  Figure  10.  The  first  step  in  the  procedure  is  to 
obtain  a  low  spatial  resolution  measurement  of  the  entire  velocity  field  using  a  constant  spatial 
shift  between  all  first  and  second  interrogation  windows  at  every  measurement  point.  This  coarse 
measurement  provides  an  initial  estimate  of  the  velocity  field.  In  the  case  of  the  Couette  flow 
measurements,  64  pixel  square  interrogation  windows  with  a  spatial  shift  between  first  and 
second  interrogation  windows  equal  to  the  estimated  mean  displacement  of  the  particles  in  the 
images  was  used.  No  overlap  between  neighboring  windows  was  used.  After  the  velocity  field 
was  calculated,  the  vectors  were  validated  using  CleanVec1.  CleanVec  compares  a  candidate 
measurement  with  the  median  of  its  neighbors’  values  to  determine  whether  it  is  more  likely  a 
valid  measurement  or  an  invalid  measurement  attributable  to  noise.  If  a  vector  differs  from  the 
median  of  its  neighbors  by  more  than  an  amount  specified  by  the  user,  it  will  be  replaced  by 
either  the  second  choice  or  third  choice  correlation  peak  which  are  found  along  with  the  primary 
correlation  peak  and  stored  for  error  correction  purposes.  If  neither  of  these  secondary  peaks 
matches  the  neighborhood  any  better,  a  vector  is  calculated  by  interpolating  from  the  neighbors’ 
values.  Gaussian  smoothing  is  used  to  remove  noise  in  the  displacement  field. 

Using  the  initial  validated  velocity  field  measurement,  a  grid  of  new  measurement  points  is 
chosen  for  subsequent  high  spatial  resolution  interrogation  of  the  images.  A  bicubic 
interpolation  scheme  is  used  to  resample  the  initial  low  resolution  velocity  field  at  each  new  high 
resolution  measurement  point.  This  step  is  used  purely  as  an  initial  estimate  of  the  high 
resolution  velocity  field  and  is  the  only  interpolation  step  in  the  algorithm.  When  using  a 
forward  difference  formula  to  approximate  the  velocity,  each  measurement  point  is  located  at  the 
centroid  of  each  first  interrogation  window.  In  the  central  difference  formulation,  however,  the 
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measurement  point  is  located  halfway  between  the  centroids  of  the  first  and  second  interrogation 
windows.  At  each  measurement  point  in  the  CDI  algorithm,  the  centroid  of  the  first 
interrogation  window  is  displaced  from  the  location  of  the  measurement  point  by  the  integer 
number  of  pixels  that  is  nearest  -0.5*  AX,  i.e.  opposite  to  the  flow  direction,  while  the  centroid  of 
the  second  interrogation  window  is  displaced  from  the  first  by  the  integer  number  of  pixels 
nearest  AX.  Here  AX  represents  the  local  expected  displacement  measurement  from  a  previous 
iteration.  The  grid  of  measurement  points  is  chosen  so  that  the  maximum  area  of  the  images  is 
interrogated  without  any  of  the  interrogation  windows  going  beyond  the  edges  of  the  images. 

The  new  high  resolution  displacement  field  is  then  calculated  using  the  interpolated  displacement 
field  to  determine  the  spatial  shifting  of  the  interrogation  windows.  After  each  iteration,  the 
vectors  are  checked  for  convergence  to  see  if  the  velocity  field  is  progressing  to  a  value  that  does 
not  change  between  iterations.  If  another  iteration  is  needed,  the  data  are  validated  again.  If  not, 
the  data  are  usually  of  high  enough  quality  that  further  post-processing  is  not  necessary. 

Application  to  Cylindrical  Couette  Flow 

In  order  to  study  the  relative  accuracy  of  different  PIV  algorithms,  it  is  necessary  to  compare  the 
results  from  interrogating  a  known  velocity  field  using  each  of  the  different  schemes.  The 
scheme  producing  results  that  most  closely  match  the  known  velocity  field  will  be  the  more 
accurate  one.  Two  adaptive  image  shifting  techniques  will  be  compared.  The  first  is  a 
conventional  forward  difference  interrogation  (FDI)  scheme  and  the  second  is  the  new  central 
difference  interrogation  (CDI)  scheme.  These  two  adaptive  techniques  both  require  only  two 
images  to  compute  and  approximately  the  same  amount  of  computational  effort  when  compared 
to  each  other  but  do  require  additional  computational  effort  when  compared  with  constant  shift 
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techniques.  Cylindrical  Couette  flow  is  used  as  a  test  bed  to  compare  the  two  PIV  algorithms 
because  it  is  a  relatively  simple  flow  with  a  closed-form  solution  in  which  experimental  velocity 
measurements  can  be  compared  to  Monte  Carlo  simulated  images  and  analytical  calculations  to 
determine  the  relative  accuracy  of  each  technique. 

Cylindrical  Couette  tlow  is  produced  between  two  differentially  rotating  concentric  cylinders. 

The  steady  state  velocity  field  is  nearly  one-dimensional  far  from  the  ends  of  the  apparatus  and 
obeys  the  simple  and  well  known  solution  to  the  Navier-Stokes  equations 

vr=vz=0  (1: 

ve(r)=r-w(r)=  A0  -r  +  B0/r 

where  vr,  v9,  and  vz  are  the  radial,  azimuthal,  and  axial  velocity  components,  respectively,  r  is  the 
radial  position,  ©  is  the  angular  velocity,  and  the  no-slip  boundary  conditions  at  the  inner  and 
outer  cylinder  uniquely  determine  the  constants  A0  and  B0.  Since  the  velocity  field  is  known,  a 
flow  can  be  simulated  by  evolving  randomly-distributed  particles  forward  in  time  using  the 
velocity  field  in  Equation  (15).  Wereley  and  Lueptow  (1994,  1998),  among  others,  have  shown 
that  experimental  me'Wdrements  of  this  flow  very  closely  approach  the  analytical  solution, 
making  it  a  good  candidate  for  comparing  simulated  flows  to  experimental  results. 

Forward  Difference  Error  Evaluation 

Because  PIV  is  usually  performed  on  a  Cartesian  grid  while  cylindrical  Couette  flow  exhibits 
azimuthal  symmetry,  notation  can  become  somewhat  confusing.  To  bridge  this  notational  gap, 
we  will  consider  a  fluid  element  located  at  the  Cartesian  coordinates  (X0,Yo)  which  coincide  with 
the  polar  coordinates  (r0,  Go)  (see  Figure  1  la).  The  angular  velocity  of  the  fluid  element  at  a 
position  r0  is  given  by  ©(ro).  This  fluid  element  will  be  swept  through  an  arc  of  length  r0(o(ro)At 


during  the  time  At  to  the  new  position  (Xi,Y i)  which  coincides  with  (ro,0i).  Because  the  velocity 
field  is  independent  of  azimuthal  position,  we  can,  without  loss  of  generality,  set  90=0  for  the 
remainder  of  the  analysis  to  simplify  the  mathematics.  The  standard  first  order  accurate  forward 
difference  discretization  of  the  velocity  field  is 

V  X(t  + At)— X(t)  n 

’  meas  ^  ' 

Since  the  fluid  elements  trace  out  curved  paths,  and  this  discretization  is  a  two  point,  linear 
approximation  to  the  derivative,  there  will  necessarily  be  some  discrepancy  between  the 
measurement  and  the  actual  velocity.  In  order  to  compare  the  two  PIV  algorithms,  we  will 
quantify  this  error  for  both  schemes.  The  Cartesian  components  of  the  velocity  field  can  be 
written  as 


Xx-X0  r0  cos  {(O  (r0  )At}-r0 


_  ^0  __ 

V  x,meas  ^ 


At 


Yx-Y0  /-0sin{ft)(r0)At} 


l\ _ to 

At 


At 


(17a) 

(17b) 


The  actual  instantaneous  velocity  at  the  point  f'Xj,Y0)  can  be  written  as 

V  =0 

vx,act  u 
^7y,act  —  rO®(ro) 


(18) 


The  difference  between  these  two  velocities  will  tell  us  how  much  of  an  effect  the  limitations  of 
the  PIV  algorithm  will  have  on  the  overall  accuracy  of  the  measurements.  To  simplify  the 
resulting  expression,  the  trigonometric  functions  will  be  expanded  as  power  series.  After 
significant  algebraic  manipulation,  the  error  in  the  X  and  Y  direction  are  found  to  be 
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(19) 


AVX  =  -  ^  r0  At  ■  co(r0  )2  +  ^  r0  At 3  co(r0  )4  + . . . 

AVy  =-gr0At2O)(r0)3  +^r0At4co(ro)5  +••• 

where  AVX  and  AVy  represent  the  bias  error  in  each  direction. 

The  error  in  X  (or  radial)  component  is  proportional  to  At  while  the  error  in  the  Y  (or  azimuthal) 
component  is  proportional  to  At2.  In  a  general  flow  with  no  particular  symmetry,  the  error  in 
both  components  should  be  order  At,  but  the  azimuthal  symmetry  in  this  flow  increases  the 
accuracy  of  the  azimuthal  component  to  order  At2.  In  the  special  case  of  measuring  Couette 
flow,  it  is  not  necessary  to  measure  the  radial  velocity  component  since  it  should  be  zero  (unless 
to  verify  that  the  experiment  is  properly  aligned).  However,  in  cases  where  higher-order  flow 
states  are  being  measured,  such  as  Taylor  vortex  flow  and  wavy  Taylor  vortex  flow,  where  the 
radial  velocity  component  is  of  interest  and  is  much  smaller  than  the  azimuthal  component,  first 
order  accuracy  in  At  can  be  restrictive. 

Central  Difference  Error  Evaluation 

Although  the  CDI  technique  is  ostensibly  second-order  accurate  in  At,  the  accuracy  of  the 
technique  must  be  carefully  examined  to  prove  its  advantages.  Toward  this  end,  we  will 
consider  a  fluid  element,  initially  at  (X.|/2,Y.i/2),  corresponding  to  (r0, 6-1/2).  which  moves  with  the 
fluid  surrounding  it  to  (X+i/2,Y+i/2),  corresponding  to  (ro,0+i/2)  (see  Figure  1  lb).  Consequently, 
the  measurement  point  at  (Xo,Yo)  is  no  longer  on  the  path  followed  by  the  fluid  element.  The 
components  of  the  velocity  field  can  be  written  as 

v,.„  =  ^  =  0,d  (20) 
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(21) 


V 


y,meas 


Y+l/2  -Y-l/2 
At 


At 


The  actual  instantaneous  velocity  at  the  point  (Xo,Yo)  can  be  written  as 


V  =0 

v  x,act  u 


Vy,act  =  ro  cos|(o(r0 ) y  jw  r0  COS  jco(r0  )y 


(22) 


The  difference  between  the  actual  velocity  and  measured  velocity  can  be  better  assessed  after 
expanding  both  solutions  as  a  power  series  and  combining  terms  of  the  same  order  in  At 


where  AVX  and  AVy  represent  the  bias  error  in  each  direction. 

Using  the  central  difference  formulation  for  the  velocity  approximation,  we  would  expect  the 
error  in  both  velocity  components  to  be  on  the  order  of  At2.  However,  there  is  no  error  in  the 
radial  velocity  component  due  to  the  azimuthal  symmetry  of  the  flow.  The  error  in  the  azimuthal 
component  is  order  At2  as  expected. 


Monte  Carlo  Simulations 

One  method  for  evaluating  the  accuracy  of  PIV  algorithms  is  to  process  a  set  of  simulated 
particle  images  displaced  according  to  a  prescribed  velocity  field.  This  approach  has  certain 
advantages  over  evaluating  algorithms  using  experimental  results,  which  will  naturally  have 
errors  and  inaccuracies  associated  with  them.  Sets  of  simulated  images  can  easily  be  constructed 
for  a  variety  of  operating  parameters  by  changing  a  few  programming  parameters.  This 
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approach  was  used  to  assess  the  importance  of  operating  parameters  in  the  PIV  process  by 
Guezennec  and  Kiritsis  (1990)  and  Keane  and  Adrian  (1990)  and  many  subsequent  studies. 
Image  pairs  are  simulated  by  creating  a  randomly  distributed  particle  image  field  (hence  the  use 
of  the  term  ‘Monte  Carlo’)  in  the  flow  domain  and  then  evolving  the  positions  of  those  particles 
forward  in  time  according  to  the  prescribed  velocity  field;  in  this  case  Equation  (15)  above. 
Based  on  the  observation  of  many  real  particle  images,  the  particle  image  intensity  can  be 
modeled  as  a  Gaussian  curve.  Because  CCD  cameras  are  composed  of  many  small  pixels,  the 
Gaussian  distribution  of  image  intensities  should  be  averaged  over  the  area  of  each  pixel  as 
demonstrated  by  Huang,  et  al.  (1997).  While  this  step  is  very  important  when  investigating  the 
effect  of  subpixel  resolution  interpolation  routines,  it  is  not  important  in  the  current  investigation 
because  any  imprecision  in  the  particle  image  shape  will  be  constant  among  the  algorithms  being 
evaluated. 

In  order  to  facilitate  direct  comparisons  between  the  simulations  and  the  experiments,  the 
parameters  of  the  simulations  were  closely  matched  to  the  experimental  parameters.  Even 
though  the  important  parameters  were  matched,  the  simulated  images  were  idealized  in  that  all 
the  particle  images  had  the  same  diameter  and  the  same  intensity.  All  of  the  simulated  particle 
images  were  set  to  the  same  maximal  gray  level  intensity  (107)  that  matched  the  mean  particle 
brightness  above  the  background  level  in  the  experiments.  The  background  level  of  the 
simulations  was  set  to  be  zero  since  the  background  only  contributes  to  a  DC  component  in  the 
correlation  plane.  The  mean  particle  image  diameter  (3.1  pixels)  and  the  mean  particle  image 
concentration  (3.3  particle  images  in  a  16x16  pixel  interrogation  window)  were  also  matched  to 
typical  experimental  conditions.  Twenty-six  sets  of  image  pairs  were  generated  for  each  set  of 
operating  conditions  explored. 
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Experimental  Measurements 


The  experimental  flow  system,  shown  in  Figure  12,  consisted  of  two  concentric  acrylic  cylinders, 
a  rotating  inner  cylinder  and  a  stationary  outer  cylinder  with  nominal  radii  ri=4.24  cm  and 
r0=5.23  cm,  respectively.  The  resulting  gap  width  was  d=ro-r,=0.99  cm,  and  the  radius  ratio  was 
T|=  rj/ro=0.8 1 .  The  inner  cylinder  was  driven  by  a  stepper  motor  at  a  rate  of  0. 1 06  Hz,  below  the 
transition  to  Taylor  vortex  flow.  The  stepper  motor  allowed  the  rotational  speed  to  be  precisely 
controlled.  The  actual  rotation  rate  varied  by  ±  0.001  Hz  which  was  monitored  by  an  optical 
encoder  having  a  resolution  of  300  pulses  per  revolution.  The  two  cylinders  were  held  concentric 
by  precision  bearings  mounted  in  opposing  aluminum  endcaps.  The  ratio  of  the  length  of  the 
annulus  to  the  gap  width  was  r=47.8.  The  particles  were  illuminated  by  a  laser  light  sheet 
formed  in  a  plane  perpendicular  to  the  common  axis  of  the  cylinders  and  imaged  through  an 
acrylic  window  mounted  in  the  lower  endcap.  The  working  fluid  was  a  water-glycerine  mixture, 
40%  by  weight,  with  density  1.1  g  cm’3  and  viscosity  3.5  cSt  at  21. 5°C.  The  flow  was  seeded 
with  16  pm  silvered  hollow  glass  spheres  having  a  density  of  1.6  g  cm'3.  The  particles  remained 
suspended  due  to  their  small  size  and  have  been  shown  to  accurately  follow  the  flow  (Wereley 
and  Lueptow,  1998). 

Two  25  mJ  per  pulse  Nd:YAG  lasers  operating  at  20  Hz  were  used  to  illuminate  the  flow  at 
consecutive  time  instants.  The  laser  pulses  were  separated  by  7.5,  15,  30,  and  60  milliseconds, 
depending  on  the  run,  to  study  the  dependence  of  the  measurement  error  on  the  time  step.  The 
laser  sheet  had  a  thickness  of  approximately  1  mm,  with  ±0.5%  variation  in  thickness  in  the 
measurement  field.  It  was  carefully  aligned  in  a  radial-azimuthal  plane  to  avoid  crosstalk 
between  velocity  components.  Pairs  of  images  each  containing  particle  images  from  a  single 
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laser  flash  were  obtained  using  a  TSI  cross-correlation  CCD  camera,  positioned  perpendicular  to 
the  laser  sheet.  A  Computar  telecentric  lens  (model  TEC-M55)  with  a  55  mm  focal  length  and  an 
aperture  of  £2.8  with  2X  extender  was  used  to  avoid  distortion  caused  by  parallax.  Twenty-six 
sets  of  image  pairs,  the  same  number  of  image  sets  generated  in  the  Monte  Carlo  investigation, 
were  recorded  for  each  set  of  operating  conditions  explored. 

Velocity  Field  Comparison 

The  sets  of  synthetic  images  produced  through  the  Monte  Carlo  technique  (discussed  in  Section 
4.2)  and  the  experimental  images  acquired  (discussed  in  Section  4.3),  were  each  analyzed  twice. 
The  first  analysis  was  the  forward  difference  interrogation,  while  the  second  was  the  central 
difference  interrogation.  While  studying  the  convergence  of  the  iterative,  adaptive  algorithms  is 
not  the  primary  purpose  of  the  present  work,  it  is  certainly  an  important  part  of  it.  After  each 
iteration  is  completed,  the  velocity  produced  from  that  iteration  is  used  as  the  adaptive  shifting 
file  for  the  next  iteration  without  alteration  or  correction.  In  general,  the  velocity  data  could  be 
validated  or  smoothed  at  this  point  but  that  was  not  necessary  for  these  measurements.  After 
repeating  this  procedure  several  times,  the  measurement  at  a  given  point  will  either  have  reached 
a  final  value  which  does  not  change  with  subsequent  iterations  or  it  will  be  oscillating  back  and 
forth  between  two  or  more  possible  solutions.  This  oscillation  is  caused  by  the  adaptive  shift  in 
one  iteration  leading  to  a  measurement  that  does  not  produce  the  same  shift  in  the  next  iteration. 
While  this  oscillation  might  at  first  seem  troubling,  it  is  very  easy  to  explain  and  quantify. 

Figure  13  shows  two  measures  of  the  convergence  rate  of  the  CDI  algorithm  as  applied  to  the 
Taylor-Couette  measurements  here.  The  first  measure  is  how  quickly  the  RMS  difference 
between  the  measured  velocity  and  actual  velocity  approaches  a  final  value  and  the  second 
measure  is  the  fraction  of  measurement  points  have  NOT  converged.  Both  of  these  measures  are 
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plotted  in  Figure  13  as  a  function  of  iteration  number.  Both  are  also  well-modeled  by  a  decaying 
exponential  which  is  superposed  on  each  set  of  data  point  symbols.  This  figure  demonstrates 
that  with  only  a  single  iteration,  the  measurements  are  accurate  to  within  1%  of  the  correct  value, 
while  after  3  iterations,  the  measurements  are  accurate  to  within  0.1%,  and  after  5  iterations,  the 
measurements  are  accurate  to  within  0.02%,  which  is  as  accurate  as  they  will  ever  get  in  these 
experiments.  Figure  13  also  shows  that  after  4  iterations,  all  but  1%  of  the  measurements  have 
reached  a  steady  state  value  after  which  further  iterations  will  produce  no  additional 
convergence.  The  iteration  scheme  used  in  this  work  is  a  very  simple  one,  merely  feeding  back 
the  velocity  field  produced  in  one  iteration  to  calculate  a  more  accurate  velocity  field  in  the  next 
iteration,  which  worked  adequately  in  this  situation.  A  more  robust  approach  would  be  to 
include  some  damping  in  each  iteration  to  prevent  the  measurements  from  running  away  from 
the  correct  solution  due  to  a  bad  measurement  in  one  iteration  and  to  prevent  the  oscillating 
measurements  mentioned  above.  In  addition,  the  fundamental  assumption  of  PIV,  that  the 
particle  images  be  randomly  distributed,  makes  running  away  from  the  correct  measurements 
statistically  unlikely.  As  a  safety  check,  any  CDI  measurement  that  differs  from  the  FDI 
measurement  made  at  the  same  point  by  more  than  some  threshold  amount  can  be  replaced  with 
the  FDI  measurement. 

Because  the  experimental  conditions  did  not  exactly  match  the  Monte-Carlo  simulation 
conditions,  the  data  are  expressed  as  dimensionless  quantities  to  facilitate  comparisons  between 
the  two  data  sets.  The  radius  ratio  r\  was  matched  between  the  experiments  and  the  simulations 
while  the  actual  size  (in  pixels)  of  the  inner  and  outer  cylinders  differed  somewhat.  For  the 
simulations,  the  outer  cylinder  radius  measured  890.0  pixels  while  for  the  experiments  it 
measured  23 1 6. 1  pixels.  Particle  images  were  similarly  larger  in  the  experiments  so  that  the 
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interrogation  windows  used  in  analyzing  the  experiments  were  square,  measuring  16  pixels  on  a 
side,  while  for  the  simulations  they  were  squares  measuring  8  pixels  on  a  side.  The  errors  are 


expressed  as  fractional  errors  relative  to  the  linear  speed  of  the  surface  of  the  inner  cylinder  ve.ic 
since  that  is  the  only  velocity  scale  available.  The  time  delay  between  images  At  was  expressed 
nondimensionally  as  the  angular  rotation  measured  in  radians  between  the  two  laser  flashes. 
Using  these  two  dimensionless  quantities,  the  radial  velocity  measurement  error  8vr  in  the 
forward  difference  formulation  expressed  in  Equation  (19)  can  be  rewritten  as 


ve,ic 


2  V9(ric) 


(24) 


where  5vr/ve>jC  is  the  fractional  error  and  to(r)  angular  rotation  rate  of  the  fluid  at  radius  r  which  is 
given  in  Equation  (15).  Q*(r)  is  the  modified  angular  displacement  between  laser  flashes, 
measured  in  radians,  of  a  fluid  element  at  a  radial  position  r.  This  modified  angular  displacement 
varies  from  zero  at  the  outer  cylinder  to  <T1C  at  the  inner  cylinder.  In  between  the  two  cylinders, 
the  modified  angular  displacement  is  scaled  by  the  ratio  of  the  local  fluid  velocity  to  the  inner 
cylinder  surface  speed.  From  this  relationship  it  is  clear  that  the  fractional  error  of  the  forward 
difference  scheme  will  depend  both  on  the  radial  position  as  well  as  on  the  time  between  laser 
flashes  At.  The  parameter  Q  can  have  several  interpretations.  It  can  be  thought  of  as  a 
dimensionless  time  step  At  because  it  is  directly  proportional  to  At  for  a  fixed  radial  position.  It 
can  also  be  thought  of  as  a  dimensionless  measure  of  radial  position  because  it  is  a  monotonic, 
albeit  nonlinear,  function  of  radial  position  when  the  time  step  At  is  fixed.  Consequently  it  is  an 
ideal  parameter  for  comparing  the  experimental  and  computational  results  because  it  will  allow 
both  the  effects  of  time  step  size  as  well  as  radial  position  to  be  explored. 
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To  study  the  dependence  of  algorithm  accuracy  on  increasing  time  delay  lengths  At,  a  succession 
of  time  delays  and  radial  positions  were  studied.  After  the  velocity  fields  were  calculated  for 
each  data  set,  the  Cartesian  velocity  vectors  were  decomposed  into  polar  velocity  vectors  based 
on  the  location  of  the  inner  and  outer  cylinder  walls.  Any  radial  velocity  component  signifies  an 
error  because  the  flow  should  be  purely  azimuthal.  The  azimuthal  error  was  found  by 
subtracting  the  measured  azimuthal  component  from  that  predicted  by  Equation  (15).  Since  the 
azimuthal  error  is  order  At2  for  both  algorithms,  it  does  not  provide  any  insight  into  which 
algorithm  will  be  more  effective  and  will  not  be  discussed  here.  However,  the  radial  error  should 
be  proportional  to  At  as  well  as  a  function  of  radial  position,  as  Equation  24  shows,  for  the 
forward  difference  algorithm  while  the  error  should  not  depend  on  the  time  delay  for  the  central 
difference  algorithm. 

The  dependence  of  the  radial  velocity  error  on  radial  position  was  explored  first.  The  radial 
velocity  component  was  divided  by  the  surface  speed  of  the  inner  cylinder  and  plotted  as  a 
function  of  both  the  parameter  £2*  (bottom  axis)  as  well  as  the  radial  position  (top  axis)  in  Figure 
14.  For  a  fixed  inner  cylinder  speed  and  the  relatively  large  radius  ratio  of  this  experimental 
apparatus  (r|=0.8 1 ),  the  parameter  £2*  is  nearly  a  linearly  function  of  radial  position  with  £2 
equal  to  zero  at  the  outer  cylinder,  a  radial  position  equal  to  1 .00,  and  £2  equal  to  the  angular 
displacement  of  the  inner  cylinder  between  laser  flashes  at  the  inner  cylinder,  a  radial  position 
equal  to  the  radius  ratio  t|=0.81.  In  Figure  14,  each  of  the  ‘+’  symbols  represents  a  single 
velocity  measurement  point  while  the  solid  line  represents  the  predicted  —1/2  slope  line.  The 
dashed  lines  represent  one  standard  deviation  (±0.0048)  about  the  -1/2  slope  line.  The 
concentration  of  the  symbols  is  greatest  at  the  left  side  of  the  figure — the  side  of  the  figure  that 
corresponds  to  the  outer  cylinder — and  smallest  at  the  right  side  of  the  figure — the  side  of  the 
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figure  that  corresponds  to  the  inner  cylinder.  Note  that  the  symbols  are  not  distributed  all  the 
way  to  the  inner  or  outer  cylinders  because  the  finite  size  of  the  interrogation  region  excludes  the 
centroids  of  the  measurement  points  from  reaching  either  wall.  The  concentration  of 
measurement  points  varies  significantly  because  of  the  area  effect,  i.e.  larger  radial  positions 
correspond  to  larger  areas  containing  more  measurement  points.  The  general  sloping  distribution 
of  the  measurement  points  confirms  the  prediction  about  how  the  error  depends  on  radial 
position  while  the  scatter  in  the  measurement  points  (±0.0048)  demonstrates  the  size  of  the 
random  error  in  relation  to  the  systematic  bias  error.  Near  the  outer  cylinder  the  random  error 
dominates  while  near  the  inner  cylinder  the  bias  error  dominates. 

To  quantify  the  bias  error,  i.e.  radial  velocity  component,  as  a  function  of  increasing  time  delays 
At,  for  each  At  the  radial  velocity  measurements  are  plotted  as  a  function  of  radial  position  (as  in 
Figure  14)  and  a  second  order  polynomial  is  fitted  to  the  measurements  to  allow  for  curvature 
effects.  The  radial  velocity  that  would  be  measured  at  the  surface  of  the  inner  cylinder  (where 
the  bias  error  is  worst)  is  calculated  by  finding  the  polynomial’s  value  at  the  inner  cylinder  radial 
position.  This  quantity  is  plotted  in  Figure  15  as  a  function  of  angular  displacement  of  the  inner 
cylinder  Q*.  Figure  15  shows  that  for  the  forward  difference  interrogations  of  both  the 
experimental  images  as  well  as  the  simulated  images,  the  fractional  errors  calculated  at  each 
cylinder  rotation  Q.*  (dimensionless  At)  fall  very  near  the  predicted  -1/2  slope  line.  Similarly, 
the  central  difference  interrogations  for  both  the  experimental  and  simulated  images  are  scattered 
about  the  predicted  error  of  zero.  We  would  expect  these  trends  to  continue  until  other  error 
terms  begin  to  dominate,  such  as  the  spatial  gradients  associated  with  larger  inner  cylinder 
rotations. 
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3.  Application  to  Micronozzles 

There  is  a  significant  effort  in  the  aerospace  community  to  miniaturize  aircraft  and 
spacecraft.  In  recent  years,  DARPA  has  focused  resources  on  developing  aircraft  that  are  less 
than  13  cm  in  length.  Miniaturized  aircraft  and  spacecraft  may  potentially  utilize  a  variety  of 
microfluidic  devices  such  as  micro-turbines  for  power  generation,  microthrusters  for  thrust 
vectoring,  and  micronozzles  for  aerodynamic  control. 

Since  the  length  scales  of  these  devices  are  typically  on  the  order  of  tens  of  microns,  the 
Reynolds  numbers  that  characterize  these  flows  are  several  orders  of  magnitude  lower  than 
comparable  devices  at  the  macroscale.  The  Mach  numbers  associated  with  these  microdevices 
are  typically  the  same  order  of  magnitude  as  their  macroscopic  counterparts.  Relatively  little 
attention  has  been  focused  on  understanding  flows  with  high  Mach  numbers  and  low  Reynolds 
numbers.  A  notable  exception  is  the  work  by  Harley,  et  al.  (1995)  and  others. 

The  next  century  of  space  exploration  will  emphasize  the  use  of  spacecraft  that  are  less 
expensive,  relatively  simple,  and  significantly  smaller  in  scale  than  current  spacecraft  (Tang, 
1998).  Advancements  in  microsystems  technology  will  certainly  play  a  role  in  the  development 
of  these  new  spacecraft.  Microspacecraft  require  high  precision  thrust  control  for  maneuvering 
and  attitude  control.  For  example,  the  thruster  requirements  for  a  10  Kg  spacecraft  is  to  deliver 
impulse  bits  (I-bits)  of  5  pN  s  with  thrust  levels  of  1.8  mN.  According  to  Tang  [2],  the  World’s 
smallest  rocket  engine  is  currently  the  Moog  58x  125,  which  delivers  an  I-bit  of  100  pN  s  at  the 
thrust  level  of  4.5  mN.  To  meet  these  design  criteria,  more  precise  microthrusters  must  be 
developed.  Characteristically,  cold  gas  thrusters  produce  lower  impulse  bits  than  what  can  be 
produced  by  more  efficient  thrusters  such  as  resistojets  and  arcjets  (Mueller,  1997). 
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Nearly  all  the  supersonic  nozzles  that  are  currently  used  in  spacecraft  have  been  designed 
using  inviscid  flow  calculations.  However,  at  the  microscale,  Reynolds  numbers  are  relatively 
low  and  viscosity  effects  cannot  be  neglected.  The  optimal  nozzle  shape  and  throat  size  for  a 
specific  thrust  level  remains  unclear.  Investigations  of  micronozzle  performance  have  been 
limited  to  numerical  simulations,  or  measurement  of  bulk  flow  properties  (Mueller,  1997;  Kim, 
1994;  Janson  and  Heelvajian,  1996;  Bayt  et  al.,  1997). 

Bayt  et  al.  (1997)  fabricated  micronozzles  by  Deep  Reactive  Ion  Etching  (DRIE)  2-D  nozzle 
contours,  entrance  tanks,  and  exit  tanks  into  500  pm  thick  Silicon  wafers.  An  example  of  this 
type  of  micronozzle  is  shown  in  Figure  16.  The  top  and  bottom  of  the  Silicon  wafers  are 
anodically  bonded  to  500  pm  thick  glass  wafers.  Holes  for  fluid  input  and  output  are  drilled  into 
one  of  the  glass  wafers. 

Micron  Resolution  Velocity  Measurements 

Until  recently  there  has  been  no  well-utilized  technique  for  measuring  the  details  of  flow 
fields  inside  microfluidic  devices.  Meinhart  et  al.  (199  gave  a  review  of  the  available  diagnostic 
techniques  for  micro  fluidics.  Initial  work  in  this  area  was  conducted  by  Lanzillotto  et  al.  [8], 
They  used  x-ray  micro-imaging  to  record  the  motion  of  1  -  20  pm  emulsion  droplets  mixed  in  a 
liquid  flow.  The  primary  advantage  of  the  x-ray  technique  is  that  it  does  not  require  optical 
access  to  image  inside  the  microfluidic  devices.  However,  it  is  limited  to  relatively  slow  flows 
and  has  relatively  low  spatial  resolution. 

The  first  successful  micron-resolution  PIV  experiment  was  reported  by  Santiago  et  al. 
(1998).  They  used  300  nm  dia.  fluorescing  flow-tracing  particles  to  measure  a  steady  Hele-Shaw 
flow  around  a  30  pm  wide  obstacle,  with  approximately  7  pm  x  7  pm  x  1.5  pm  spatial 
resolution.  Meinhart  et  al.  (1999a)  extended  the  technique,  using  a  pair  of  Nd:YAG  lasers  to 
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illuminate  200  nm  dia.  fluorescing  particles,  to  measure  the  flow  through  a  30  |tm  x  300  pm 
channel.  A  velocity-vector  field  measured  in  the  microchannel  is  shown  in  Figure  3.  The  spatial 
resolution  of  the  measurements  were  13.8  pm  x  0.9  pm  x  1.8  pm.  The  velocity  vectors  were 
overlapped  by  50%,  yielding  a  velocity-vector  spacing  of  450  nm  in  the  wall-normal  direction. 

The  micro  Particle  Image  Velocimetry  (PIV)  technique  presented  here  is  the  first  step 
towards  measuring  detailed  flow  fields  inside  micronozzles.  There  exist  significant  challenges 
that  must  be  overcome  in  order  to  measure  accurately  velocity  fields  of  high-speed  flows  at  the 
microscale.  These  measurements  will  provide  information  about  subsonic  boundary  layer 
growth,  flow  separation,  and  the  streamwise  evolution  of  the  supersonic  core  flow. 


Flow-Tracing  Particles 

The  desire  to  measure  high-speed  velocity  fields  accurately  inside  micronozzles  with  order 
one-micron  spatial  resolution  places  significant  constraints  on  the  choice  of  flow-tracing 
particles.  The  particles  must  be  significantly  smaller  than  the  length  scale  of  the  microdevice, 
and  must  be  small  enough  to  faithfully  follow  the  flow  under  extremely  large  accelerations. 

If  a  supersonic  flow  through  a  micronozzle  accelerates  from  Mach  1  at  the  throat  to  Mach.2 
one  hundred  microns  downstream,  and  assuming  uniform  acceleration,  the  flow  will  accelerate  at 
an  average  rate  of  approximately  9.0  E  +08  m  s'2,  or  approximately  100  million  times  the 
acceleration  of  gravity.  Assuming  Stokes’  drag,  the  difference  between  the  particle  velocity  and 
the  fluid  velocity  can  be  estimated  by 


v-u| 


36pv 


(25) 


where  v  is  the  particle  velocity,  u  is  the  velocity  of  the  fluid,  v  represents  the  total  time 
derivative  of  the  particle  velocity,  and  the  subscript  p  represents  the  properties  of  the  particle 
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(Adrian,  1991).  Using  Silica  (Si02),  which  has  a  density  of  2  g  cm'3,  a  100  nm  dia.  particle  will 
slip  ~  27  m  s'1,  while  a  50  nm  dia.  particle  will  slip  ~  7  m  s'1.  This  represents  about  4  -  8%  FS 
error  for  the  100  nm  and  1  -  2%  FS  error  for  50  nm  dia.  particles.  Under  these  conditions,  a  100 
nm  dia.  particle  represents  an  upper  bound  on  the  size  of  particle  that  should  be  used  for  PIV 
measurements  of  high-speed  flow  in  micronozzles. 

The  slip  flow  calculation  from  (25)  assumes  Stokes’  drag  is  valid,  which  is  a  result  of 
continuum  flow  theory.  The  mean-free  path  of  N2  at  10  atmospheres  is  approximately  8  nm,  so 
the  Knudsen  numbers,  Kn  =  A jdp ,  of  the  particles  at  the  throat  of  the  nozzle  are  approximately 
Kn  =  0.16  and  0.08,  for  the  50  and  100  nm  dia.  particles,  respectively.  These  Knudsen  numbers 
are  similar  to  Knudsen  numbers  of  0.5  -  1.0  pm  dia.  particles  commonly  used  in  PIV 
experiments  of  gaseous  flows  at  standard  laboratory  conditions. 

Flow-tracing  particles  must  also  be  large  enough  to  scatter  sufficient  light  so  that  their 
images  can  be  recorded.  In  the  Rayleigh  scattering  regime,  where  the  particle  diameter  is  much 
smaller  than  the  wavelength  of  light,  dp«X,  the  amount  of  light  scattered  by  a  particle  varies 

as  d~6  (Bom  and  Wolf,  1997).  Since  the  diameter  of  the  flow-tracing  particles  must  not  be  larger 

than  50  nm  -  100  nm,  in  order  to  follow  the  flow  through  supersonic  nozzles,  their  diameters  are 
1/10  to  1/5  smaller  than  the  wavelength  of  green  light,  A  =  532  nm,  and  are  therefore 
approaching  the  Rayleigh  scattering  criteria.  This  places  significant  constraints  on  the  image 
recording  optics,  making  it  extremely  difficult  to  record  particle  images. 

One  solution  to  the  imaging  problem  is  to  use  epi-fluorescence  imaging  to  record 
inelastically-scattered  light  from  fluorescently-labeled  particles  through  an  optical  filter,  which 
removes  the  background  light.  This  technique  was  used  successfully  in  liquid  flows  to  record 
images  of  200  -  300  nm  diameter  fluorescent  particles  (Santiago  et  al.,  1998;  Meinhart  et  al.. 
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1999a).  While  fluorescently-labeled  particles  are  well  suited  for  micro-PIV  studies  in  liquid 
flows,  they  are  not  applicable  to  high-speed  air  flows  for  several  reasons.  First,  commercially 
available  fluorescently-labeled  particles  are  available  only  in  aqueous  solutions.  In  principle,  the 
particle-laden  solutions  can  be  dried,  and  the  particles  emitted  into  an  air  stream.  Unfortunately, 
we  have  not  been  able  to  dry  particles  without  significant  particle  clumping.  Furthermore,  the 
emission  decay  time  of  many  fluorescent  molecules  is  on  the  order  of  several  nanoseconds, 
which  may  cause  streaking  of  the  particle  images  for  high-speed  flows.  In  contrast,  elastic 
scattering  is  effectively  instantaneous.  Dry  Silicate  (Si02)  particles  are  commercially  available 
with  diameters  ranging  between  50  -  100  nm,  and  can  be  dispersed  in  air  flows  without 
significant  clumping.  Silicate  particles  have  a  density  of  2  g  cm'3,  which  is  about  twice  the 
density  of  polystyrene  particles.  Therefore,  Silicate  particles  will  slip  twice  the  amount  of 
polystyrene  particles.  Since  the  slip  error  is  ~  d2p,  the  increase  in  error  caused  by  using  the  denser 
Silicate  particles,  compared  to  the  less-dense  polystyrene  particles,  can  be  compensated  by  using 
a  particle  that  is  V2  smaller  in  diameter. 

Spatial  and  Temporal  Resolution 

For  the  purposes  of  this  discussion,  we  will  define  temporal  resolution  as  the  duration  of 
averaging  time  required  to  obtain  a  measurement.  In  the  case  of  double  pulse  PIV,  the  temporal 
resolution  is  simply  the  time  delay  between  pulses,  At.  The  ability  to  make  accurate 
measurements  while  averaging  over  a  short  interval  in  time  is  important  in  variety  of 
experiments.  It  is  particularly  important  for  measuring  high-speed  flows  in  micronozzle 
experiments.  In  this  situation,  the  desire  to  make  PIV  measurements  on  the  micron  scale  coupled 
with  the  large  velocity  associated  with  supersonic  flows  requires  nanosecond  temporal 
resolution. 
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Consider  a  flow  near  Mach  1  at  standard  temperature,  with  a  characteristic  velocity, 
uc  =340  ms'1.  If  the  desired  particle  displacement  is,  say,  Ax  =  10/tm,  then  the  time  delay 
between  pulses  should  be,  At  ~  30  ns.  A  standard  Nd.YAG  laser  can  be  used  to  generate  a  short 
pulse  of  green  light  ( X  =  532  nm),  which  is  approximately  Gaussian  in  time  with  a  time  duration, 
8t  ~  5  ns.  During  a  5  ns  period,  a  particle  with  a  characteristic  velocity  of  uc  =  340  m  s'1  travels 
approximately  1.7  (im.  If  a  particle  moves  a  significant  fraction  of  its  effective  image  size  during 
exposure,  the  particle  image  will  be  blurred.  In  nearly  all  PIV  experiments  reported  to  date,  the 
relatively  short  light  pulse  (5t~5  ns)  from  an  Nd:YAG  laser  is  more  than  sufficient  to  freeze 
particle  motion.  It  is  only  in  the  special  case  of  high-speed  microscale  flow  that  the  pulse 
duration  of  an  Nd.  YAG  laser  becomes  limiting. 

Further  refinements  in  temporal  resolution  can  be  obtained  by  employing  mode-locked  YAG 
lasers  or  Titanium-Sapphire  (77-5)  lasers,  which  have  pulse  durations  ranging  between  50  fs  - 
100  ps.  However  for  the  current  discussion,  we  will  limit  our  investigation  to  standard  Nd.  YAG 
PIV  lasers. 

Following  the  arguments  of  MeinhaV  et  al.  [10]  the  diameter  of  the  diffraction-limited 
particle  image,  ds,  can  be  approximated  by 


d  =2.44M- 


2  NA 


(26) 


where  M  is  the  magnification,  and  NA  is  the  numerical  aperture  of  the  objective  lens. 
Approximating  the  diffraction-limited  image  and  the  geometric  image  both  as  Gaussian 
functions  one  can  estimate  the  effective  diameter  of  a  particle  image  as 

d.  =  [MV;+d;f\  (27) 
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The  blurred  image  resulting  from  a  finite  exposure  time  is  the  convolution  of  the 
instantaneous  particle  image  function  with  the  blurring  function,  where  the  blurring  function, 
B(X),  can  be  written  as 

B(X)  =  \  S(X-Xp(t))l(t)adt,  (28) 

where  X  is  the  spatial  coordinate  on  the  CCD  array,  Xv  is  the  centroid  position  of  a  particle  when 
projected  onto  the  CCD  camera,  5  is  the  Dirac  delta  function,  I  is  the  intensity  of  a  the 
illuminating  laser  light  pulse  as  a  function  of  time,  and  o  is  the  scattering  coefficient  of  the 
particles. 

Since  the  intensity  profile  of  the  laser  pulse  is  approximately  Gaussian  in  time,  the  width  of 
the  blurred  image  due  to  the  finite  exposure  time  can  be  approximated  by 

+  (»> 

The  effective  image  diameter,  de,  of  a  100  nm  dia.  particle  imaged  with  a  numerical  aperture, 
NA  =  0.6,  40x  magnification  objective  lens,  using  X  =  532  nm  light,  is  de  ~  43  pm.  If  the 
characteristic  velocity  is  uc  ~  340  m  s'1,  and  the  duration  of  the  laser  pulse  is  5t  ~  5  ns,  then  the 
width  of  the  blurred  image  is  dh ]ur  ~  80  pm.  When  projected  back  in^o  the  fluid  the  width  of  the 
blurred  image  is  ~  2  pm.  Assuming  that  one  can  measure  the  location  of  the  particle  to  within 
1/1 0th  of  width  of  the  blurred  image,  then  the  uncertainty  in  determining  the  particle  location  is  ~ 
200  nm.  For  a  relative  error  of  2%,  a  particle  image  displacement  of  approximately  10  pm  is 
required.  This  places  a  lower  limit  on  the  spatial  resolution  in  the  streamwise  direction. 

Differential  Interference  Contrast  (DIC)  Imaging 

From  the  above  discussions,  it  is  clear  that  50  -  100  nm  dia.  particles  are  required  for 
accurate  PIV  measurements  of  high-speed  flow  in  micronozzles.  Moreover,  the  time  duration  of 
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an  individual  light  pulse,  and  the  time  duration  between  light  pulses  must  be  <5t<5ns  and 
At  ~  30  ns,  respectively.  Since  the  emission  decay  time  of  most  fluorescent  molecules  is  on  the 
order  of  several  nanoseconds,  and  since  it  is  difficult  to  seed  fluorescently-dye  sub-micron 
particles  into  air  flow,  we  have  chosen  to  work  with  50  -  100  nm  dia.  dry  Silicate  particles. 

Transmitted  Differential  Interference  Contrast  (DIC)  is  a  technique  used  in  biological 
microscopy  to  image  sub-cellular  structures  in  transparent  specimens,  which  are  smaller  or 
comparable  to  the  diffraction  resolution  of  the  microscope.  Reflective  DIC  is  a  microscopic 
technique  used  in  Mineralogy  to  image  extremely  small  features  in  non-transparent  specimens. 
We  have  employed  reflective  DIC  for  the  micronozzle  experiments,  because  it  requires  optical 
access  from  only  one  direction. 

The  reflective  DIC  system  shown  in  Figure  17  consists  of  a  linear  polarizer,  a  mirror,  an 
adjustable  Wollaston  prism,  and  a  DIC  analyzer.  The  illumination  light  enters  the  system  through 
the  linear  polarizer  and  is  reflected  by  the  mirror  towards  the  Wollaston  prism.  The  Wollaston 
prism  consists  of  a  birefringent  material  oriented  at  45°  to  the  incoming  polarized  light.  The 
birefringence  causes  the  incoming  light  to  be  sheared.  The  degree  of  beam  shear  is  adjusted  by 
sliding  one  element  of  the  prism  relative  to  the  second  element.  It  is  common  to  adjust  the  shear 
so  that  it  is  less  than  the  diffraction  limit  of  the  microscope  so  that  a  double  image  is  not 
observed.  The  sheared  beam  illuminates  the  back  of  the  objective  lens,  where  it  is  refracted  into 
the  microfluidic  device.  Light  is  spectrally  reflected  from  the  surfaces  of  the  microfluidic  device, 
and  scattered  from  the  flow-tracing  particles.  The  reflected  light  from  the  particles  and  the  test 
section  is  collected  by  the  objective  lens  and  relayed  through  the  Wollaston  prism  where  the 
beam  shear  is  reversed.  Because  optical  systems  are  acausal,  light  waves  that  were  relatively 
undisturbed  in  the  test  section  will  recover  their  original  polarization,  after  reversing  back 
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through  the  Wollaston  prism.  Light  waves  that  have  been  phase  shifted  or  changed  in  amplitude 
by  objects  that  are  of  the  same  order  in  size  as  the  beam  shear  will  recover  a  different 
polarization,  after  reversing  back  through  the  Wollaston  prism.  The  light  is  then  relayed  through 
the  mirror  onto  a  DIC  analyzer.  In  the  current  configuration,  the  DIC  analyzer  is  simply  a  linear 
polarizer  oriented  90°  to  the  incoming  light.  This  filters  out  light  that  is  relatively  undisturbed  in 
the  test  section  over  length  scales  on  the  order  of  the  beam  shear.  In  practice,  the  beam  shear  is 
adjusted  to  obtain  maximum  contrast  for  a  particular  type  of  object  (Inoue  and  Spring,  1997). 

The  collected  light  is  imaged  onto  a  Princeton  Instruments  1300  x  1030  x  12  bit  cooled  CCD 
camera.  A  cooled  CCD  array  is  necessary  to  reduce  electronic  noise  to  a  sufficiently  low  level  to 
record  light  scattered  from  sub-micron  particles.  This  particular  camera  was  chosen  because  it 
has  an  interline-transfer  array  that  can  capture  two  back-to-back  images  within  a  500  ns  time 
duration.  The  time  duration  between  pulses  for  the  micronozzle  experiments  must  be  less  than 
100  nm.  Therefore,  both  pulses  of  light  must  be  recorded  on  the  same  camera  frame. 

Supersonic  flow  requires  that  the  two  light  pulses  be  separated  by  At  =  30  ns.  The  most 
common  technique  in  PIV  for  producing  two  time-sequenced  pulses  is  to  use  two  separate  lasers, 
which  are  synchronized  by  a  high-precision  electronic  pulse  generator.  Pulse  generators  are 
commercially  available  with  picosecond  resolution.  However,  Nd:YAG  lasers  typically  have  a 
jitter  of  at  least  1  ns. 

An  alternative  method  is  proposed  in  Figure  17.  It  uses  an  optical  delay  to  produce  two  time- 
sequenced  pulses  of  light.  The  primary  advantage  of  this  technique  is  that  the  pulse  sequence  is 
highly  repeatable,  and  not  subject  to  electronic  jitter.  The  primary  disadvantage  of  this  technique 
is  that  it  is  not  easy  to  change  the  time  delay  and  that  the  time  delay  cannot  be  set  to  the  precision 
of  digital  electronics. 
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The  proposed  optical  delay  consists  of  a  7J2  waveplate  that  rotates  the  vertically  polarized 
light  emanating  from  the  Nd:YAG  laser  to  45°  angle.  A  polarizing  beam  splitter  reflects  the 
vertically  polarized  component  of  the  light,  and  passes  the  horizontal  component.  Since  the 
vertical  component  has  a  longer  optical  path,  when  the  beams  are  recombine  at  the  second  beam 
splitter,  the  vertical  component  will  be  delayed  relative  to  the  horizontal  component.  A  second 
X/2  waveplate  rotates  the  horizontally  polarized  and  vertically  polarized  light  pulses  by  45°. 
Then  a  linear  polarizer  vertically  polarizes  both  light  pulses,  so  that  they  will  be  compatible  with 
the  DIC  optics. 

The  ability  of  the  DIC  system  to  image  50  nm  dia.  Silicate  particles  using  elastic  scattering  is 
demonstrated  in  Figure  18.  A  3-D  shadowing  feature  of  the  particle  images  is  characteristic  of 
DIC,  and  can  be  observed  in  Figure  18.  Figure  18a  shows  images  of  50  nm  Silicate  particles 
using  a  NA  =  1.4  oil  immersion  lOOx  objective  lens,  while  Figure  18b  shows  images  of  50  nm 
Silicate  particles  using  a  NA  =  0.6  air  immersion  40x  objective  lens.  The  particles  were  not 
observable  without  the  aid  of  DIC. 

Since  the  50  nm  particles  are  much  smaller  than  the  diffraction  resolution  of  the  microscope, 
it  is  not  clear  initially  that  individual  particles  are  being  observed  as  opposed  to  clusters  of 
particles.  To  answer  this  question,  we  took  50  nm  fluorescently-dyed  polystyrene  particles. 
These  particles  were  imaged  using  epi-fluorescence  illumination.  The  same  set  of  particles  was 
then  imaged  using  DIC.  The  epi-fluorescent  technique  utilizes  purely  in-elastic  scattering  while 
the  DIC  technique  utilizes  purely  elastic  scattering.  Therefore,  if  a  cluster  of  particles  is  imaged, 
the  particle  intensity  resulting  from  epi-fluorescent  imaging  is  proportional  to  the  number  of  dye 
molecules  being  exposed,  i.e.  the  number  of  particles  in  the  cluster.  On  the  other  hand,  the 
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particle  intensity  resulting  from  DIC  imaging  is  related  to  the  scattering  coefficient  of  the  particle 
cluster,  which  is  a  complicated  function  of  the  particle  cluster  geometry. 

The  intensity  of  the  epi-fluorescent  images  has  characteristic  quantum  levels,  corresponding 
to  a  discrete  number  of  particles  in  a  specific  cluster.  Assuming  the  smallest  increment  between 
quantum  levels  corresponds  to  the  fluorescence  of  a  single  particle,  one  can  identify  which 
clusters  contain  only  a  single  particle.  After  single  particle  images  are  identified  from  the  epi- 
fluorescence  image,  individual  particle  positions  can  be  registered  and  images  of  single  particles 
can  be  identified  on  the  DIC  image.  This  experiment  demonstrates  that  we  are  imaging  single  50 
nm  dia.  particles  using  DIC  (Matsumoto,  1999). 

Velocity  Measurements  in  a  Micronozzle 

The  micronozzles  used  in  the  water- flow  experiment  were  fabricated  by  Robert  Bayt  at  MIT. 
The  2-D  nozzle  contours  were  etched  using  Deep  Reactive  Ion  Etching  (DRIE)  in  300  pm  thick 
Silicon  wafers.  The  top  and  bottom  of  the  wafers  were  anodically  bonded  to  500  pm  thick  glass 
wafers.  Holes  were  then  drilled  into  the  glass  to  provide  input  and  output  fluid  ports.  The  wafers 
are  mounted  to  an  aluminum  manifold  and  sealed  using  O-rings  and  connected  with  plastic 
tubing  to  a  Harvard  Apparatus  syringe  pump. 

The  experiments  were  conducted  using  water  as  the  fluid  medium  to  determine  the 
conditions  in  which  PIV  measurements  could  be  obtained  in  the  relatively  deep  micronozzles. 
The  liquid  flow  was  seeded  with  relatively  large  700  nm  dia.  fluorescently-labeled  polystyrene 
particles.  The  particles  were  imaged  using  a  NA  =  0.6,  40x  objective  lens,  and  an  epi-fluorescent 
imaging  system  similar  to  that  described  by  Meinhart  et  al.  (1999a).  A  flow  rate  of  4  ml  hr  1  was 
delivered  to  the  nozzle  by  the  syringe  pump. 
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Figure  19  is  an  example  of  an  instantaneous  velocity  field  measured  in  the  throat  and 
expansion  section  of  a  straight  diffuser  with  a  24°  half  angle.  The  velocity  fields  were  calculated 
using  a  standard  cross  correlation  algorithm.  Both  interrogation  windows  were  64  x  20  pixels  in 
the  x  and  y  directions,  respectively.  When  projecting  into  the  fluid,  the  correlation  windows  were 
10.9  x  3.4  pm  in  the  x  and  y  directions,  respectively.  The  correlation  windows  were  offset  to  the 
local  particle-image  displacement  at  each  measurement  point.  The  interrogation  spots  were 
overlapped  by  50%,  yielding  a  velocity-vector  spacing  of  5.44  pm  in  the  streamwise  direction 
and  1 .7  pm  in  the  spanwise  direction. 

The  Reynolds  number,  based  upon  bulk  velocity  and  throat  width,  is  Re  =  10,  suggesting  that 
nonlinear  terms  in  the  Navier-Stokes  equation  may  be  important  in  the  expansion  region  of  the 
nozzle.  The  PIV  measurements  reveal  that  the  velocity  profile  is  nearly  parabolic  in  the  throat. 
An  inflection  in  the  velocity  profile  is  formed  as  the  flow  enters  the  expansion  section.  Further 
downstream  in  the  expansion  section,  the  flow  separates  from  the  lower  nozzle  wall.  This 
suggests  that  separation  may  also  occur  for  high-speed  gas  flows,  where  the  Reynolds  number 
can  be  two  orders  magnitude  larger  than  the  current  experiment. 


4.  Extension  to  Air  Flow 

While  the  PIV  technique  has  been  demonstrated  successfully  by  measuring  liquid  flows,  the 
measurement  of  high-speed  air  flows  through  micronozzles  remains  difficult,  in  part  because  we 
would  like  to  measure  gas  velocities  and  accelerations  much  higher  than  those  measured  in  liquid 
flows.  This  translates  into  a  requirement  for  particles  of  small  size  and  mass,  such  that  can 
faithfully  follow  the  flow.  Higher  resolution  in  time  will  also  be  necessary  to  capture  details  of 
high-speed  flow,  as  will  shorter  exposure  times  to  prevent  streaking.  These  requirements  were 
discussed  and  quantified  in  Meinhart  (1999c). 

Imaging  Solution:  DIC 

The  optical  method  of  Differential  Interference  Contrast  (DIC)  isproposed  as  a  solution  to  the 
problem  of  imaging  very  small  (50  -  100  nm)  particles  necessary  for  airflow  PIV  in 
micronozzles. 

While  DIC  has  been  found  to  be  quite  successful  in  imaging  small  particles  using  white  light, 
the  use  of  DIC  with  laser  illumination  has  not  yet  been  optimized.  Because  of  the  coherence  of 
the  laser  beam,  a  “laser  speckle”  diffraction  pattern  is  produced  when  the  DIC  system  :c.  used 
with  the  Nd:YAG  laser.  There  are  various  methods  used  to  reduce  laser  speckle.  Examples  of 
these  include:  counter-rotating  glass  diffusers  (Partlo,  and  Oldham,  1991),  modified  fiber  array 
(Dingel,  1993),  and  illumination  with  a  fluorescent  beam  produced  by  a  laser  dye. 

Seeding  the  Gas  Flow 

In  addition  to  the  restrictions  placed  on  the  airflow  PIV  system  by  the  higher  speeds  and 
accelerations,  there  is  also  the  problem  of  seeding  the  flow,  that  is,  developing  and  maintaining  a 
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sufficiently  disperse  aerosol  of  solid  or  liquid  particles  in  the  gas  stream.  There  are  at  least  four 
commonly  used  methods  for  suspending  such  an  aerosol. 

A  powder  disperser  uses  a  gas  stream,  either  aimed  at  a  powder  sample  or  bubbled  up 
through  the  powder  to  create  a  powder  cloud,  or  blown  through  a  venturi  to  suck  powder  up  into 
the  throat.  The  particles  dispersed  in  the  gas  stream  are  normally  clumps  of  several  individual 
particles  (Melling,  1997);  indeed  it  is  very  difficult  to  distribute  a  completely  deagglomerated 
aerosol.  Reducing  moisture  in  the  gas  stream  and  in  the  powder  sample  may  help  in  dispersion 
efficiency  by  reducing  hydrostatic  coagulation  forces,  but  this  must  be  balanced  with  the  effect 
of  electrostatic  forces,  which  increase  with  dryness  (Melling).  Particle  size  and  material  also 
play  a  part  in  agglomeration.  TSI  (St.Paul,  MN)  quotes  their  powder  dispersers  as  able  to 
disperse  powders  larger  than  500  nm;  anything  smaller  is  not  effectively  deagglomerated. 
Rough,  soft  particles  are  more  difficult  to  disperse  than  smooth,  hard  particles  (Melling,  1997). 

An  atomizer  avoids  the  problem  of  deagglomerating  cohesive  powders.  A  pressure 
difference  across  a  nozzle  pulls  liquid  through  the  nozzle,  creating  a  polydisperse  liquid  aerosol. 
If  the  liquid  is  a  suspension  of  (monodisperse)  solid  particles,  the  aerosol  can  be  dried,  leaving  a 
(monodisperse)  dilute  aerosol  of  solid  particles.  Alternatively,  if  the  liquid  is  a  solution  (saline  is 
common),  the  aerosol  can  be  dried,  leaving  a  solid  aerosol  of  solute  crystals. 

A  condensation  aerosol  generator  allows  a  vaporized  medium  to  condense  on  a  tiny  nucleus. 
If  this  is  done  at  constant  temperature,  as  in  a  Sinclair-LaMer  type  generator,  a  monodisperse 
aerosol  is  attained,  though  it  is  relatively  low  in  concentration. 

An  electrospray  aerosol  generator,  such  as  the  one  developed  by  TSI,  Inc.,  uses  an  electric 
field  to  pull  uniform  tiny  liquid  droplets  off  the  end  of  a  capillary  tube.  The  resulting  aerosol  is 
monodisperse  and  nearly  comparable  in  concentration  to  that  of  the  atomizer,  the  droplet  size  is 
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150  nm.  If  these  droplets  are  a  saline  solution,  final  solid  particle  size  can  be  as  large  as  50-100 
nm.  While  this  coincides  with  the  ideal  size  for  airflow  seeding  particles,  it  is  necessary  to  be 
able  to  produce  larger  (easily  seen)  particles  as  well  during  the  optimization  of  the  imaging  and 
flow  systems. 

Coagulation 

As  soon  the  aerosol  is  generated,  it  unfortunately  begins  to  recoagulate.  As  particles  move 
relative  to  each  other  by  Brownian  motion,  turbulence,  or  body  forces,  they  collide  and  coalesce. 
Theoretical  treatments  of  coagulation  begin  by  assuming  perfect  collision  efficiency,  that  is  any 
particles  which  come  into  contact  will  stick  (Green,  1964,  Fuchs,  1964,  and  Reist,  1993  all 
present  details  on  coagulation.)  If  particles  that  collide  do  not  stick  together,  they  will,  because 
of  proximity,  have  an  increased  probability  of  another  collision,  thus  raising  the  “effective” 
collision  efficiency.  If  a  monodisperse  aerosol  is  moving  only  through  Brownian  motion,  the 
number  concentration  n  as  a  function  of  time  may  be  described  as 

1  1  K  (30) 

—  - - h  Kl 

n  n0 

where  K  is  the  coagulation  constant.  K  can  be  derived  theoretically  for  the  above  case,  and  turns 
out  to  be  4,0  X  10'10  cm3/s  for  200  nm  particles  in  air  at  STP  (Fuchs,  1964),  but  additional 
complications  of  our  system  such  as  polydispersity  and  turbulence,  may  affect  this  number  by 
orders  of  magnitude.  Nevertheless,  some  interesting  observations  may  be  made  about  this 
process.  First,  for  a  high  coagulation  constant,  and  an  observation  point  a  finite  time  after  t  =  0 
(tenths  of  seconds  perhaps),  increasing  the  initial  concentration  is  meaningless:  the  particles  will 
immediately  coagulate.  There  is  therefore,  a  limit  to  the  concentration  of  an  aerosol  that  can  be 
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generated  and  maintained.  Second,  for  high  initial  concentrations,  it  is  important  for  the 
observation  point  to  be  as  soon  as  possible  after  t  =  0  so  as  not  to  allow  the  aerosol  to  coagulate. 
This  translates  for  us  into  a  short  as  possible  run  between  the  aerosol  generator  and  the  test 
section.  We  keep  these  practical  limitations  in  mind  as  we  search  for  a  flow  seeding  method. 

Airflow  Velocimetry  Experiments 

Although  it  is  common  to  use  aerosol  generators  of  the  aforementioned  type  for  seeding  gas 
flow  for  PIV,  most  experiments  have  been  done  on  much  larger  scale.  A  reduction  in  scale  not 
only  necessitates  a  reduction  in  particle  size,  but  also  an  increase  in  particle  number 
concentration.  Our  needs,  therefore,  push  the  limits  of  what  has  been  used,  what  is  available, 
and  perhaps  what  is  possible.  If  a  40X,  N.A.  0.6  objective  lens  is  used,  depth  of  field,  taken  from 
Meinhart  et  al.  (2000a),  is  10  pm,  and  resolution  of  the  camera  is  taken  to  be  6.8  pm  /  pixel,  the 
image  volume  is  3.74  X  10'7  cm3.  If  15  particles  are  needed  per  interrogation  volume  and 
interrogation  windows  of  64  X  20  pixels  are  used  as  in  Meinhart  (1999c)  corresponding  to  a 
resolution  of  5.44  X  1.7  pm,  then  it  is  desirable  to  have  about  4.5  X  105  particles  per  image 
volume.  Ensemble  averaging  technique*  for  steady  flow  can  reduce  this  nearly  an  order  of 
magnitude  (Meinhart  et  al.,  1999b),  but  that  still  leaves  a  desirable  aerosol  concentration  of  2  X 
1010  particles  /  cm3.  Specified  maximum  number  concentrations  on  commercially  available 
aerosol  generators  run  in  the  106  to  108  range,  so  dispersion  of  a  highly  concentrated  aerosol 
becomes  one  of  the  limiting  system  criteria. 

Shallow  Channel  Design 

In  order  to  ignore  for  the  moment  the  problems  associated  with  high  speed  flow,  and  to 
increase  the  overall  volume  flow  rate,  which  generally  increases  aerosol  generation,  a  new 
channel  was  designed  for  testing  the  airflow  PIV  system,  shown  in  Figure  20.  The  channel  is 
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wide,  increasing  the  flow  rate,  and  shallow,  keeping  the  illuminated  depth  to  a  minimum  to 
reduce  noise  from  out-of-focus  particles.  Additionally,  because  the  channel  is  quickly  made  on  a 
mill  rather  than  using  micromachining  techniques,  there  is  more  flexibility  than  with  the 
micronozzle.  For  example,  the  bottom  of  the  channel  was  colored  black  to  improve  image 
quality.  In  another  design,  the  channel  and  manifold  were  made  from  acrylic,  to  allow 
illumination  from  behind  the  channel.  The  intent  is  to  use  the  less  restrictive  shallow  channel  to 
get  a  PIV  system  working.  The  system  can  then  be  fine-tuned  for  the  more  restrictive 
micronozzle. 


channel  depth:  200  gm 


Figure  20.  Shallow  Channel  Design 


Powder  Dispersion  Seeder 

The  first  aerosol  generation  system  used  was  a  powder  dispersion  seeder  designed  and  built 
in-house,  diagrammed  in  Figure  21.  Nitrogen  is  the  working  gas;  TiC>2  was  chosen  as  the  powder 
because  of  its  availability  in  a  range  of  appropriate  sizes  and  its  low  density  (compared  with 
metals).  A  gas  stream  impinges  on  a  pile  of  powder,  creating  a  cloud  of  powder.  The  pressure 


differential  between  the  cloud  chamber  and  the  outside  causes  the  cloud  to  be  sucked  through  the 
channel. 

This  seeder  was  not  effective  in  producing  observable  particles  in  the  flow  stream.  Because 
there  was  clearly  powder  exiting  the  test  section,  the  absence  of  particles  in  any  images  could  be 
traced  to  either  of  two  problems:  first,  the  imaging  system  was  not  capable  of  detecting  such 
small  particles;  or  second,  the  seeder  was  distributing  such  a  low  concentration  (and  given  the 
observed  powder  mass  flow,  large  particle  size)  of  particles,  the  chances  of  one  falling  in  the 
exposure  window  were  quite  low.  Because  various  sizes  of  Ti02  particles  were  used  (200  -  700 
nm  average  diameter),  which  were  easily  imaged  when  captured  on  a  glass  slide,  and  because 
finally  one  large  (~5  (im)  particle  entrained  in  the  flow  was  captured  in  an  image,  it  seems  more 
logical  that  the  problem  is  with  the  agglomeration  of  the  powder  and  therefore  the  low  number 
concentration  of  the  aerosol.  There  was  also  obvious  loss  of  aerosol  to  tubing  walls,  which  were 
coated  with  powder  at  the  end  of  each  experiment.  Special  conductive  flexible  tubing  was  tried, 
with  no  measured  improvement  in  results. 

Atomizer 

A  TSI,  Inc.  Atomizer  (Model  3076)  was  used  in  the  airflow  velocimetry  system.  Using 
simply  water  as  a  working  liquid,  an  aerosol  of  water  droplets  was  produced.  The  shallow 
channel  quickly  became  waterlogged,  that  is,  enough  water  quickly  deposited  on  the  viewing 
window  so  as  to  prevent  imaging.  Therefore,  for  evaluation  of  the  atomizer  as  an  aerosol 
generator,  the  channel  was  removed,  and  the  flow  was  simply  aimed  at  a  glass  slide  through 
which  the  flow  was  observed.  In  this  configuration,  a  droplet  was  captured  in  an  image  about 
every  third  frame.  From  this,  an  aerosol  concentration  can  be  inferred  of  around  105  particles/cc 
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of  particles  larger  than  the  smallest  detectable  size.  TSI’s  specifications  for  this  instrument  quote 
108  particles/cc  with  a  mean  diameter  of  3.0  |im  and  a  geometric  standard  deviation  of  less  than 
2.0.  From  the  image  pairs  gathered  of  the  droplets  in  external  flow,  a  uniform  velocity  of  1.2 
m/s  was  calculated. 
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Cooled  CCD  Camera 
1300  x  1030  x  12  bit 


Figure  1.  Schematic  of  a  micro  PIV  system.  A  pulsed  Nd:YAG  laser  is  used  to 
illuminate  200  nm  diameter  fluorescent  flow-tracing  particles  through  an  epi- 
fluorescent  inverted  microscope.  A  cooled  1300  x  1030  pixel  x  12  bit 
interline-transfer  CCD  camera  is  used  to  record  the  particle  images. 


Figure  2.  Detail  of  the  nominally  30  pm  x  300  pm  rectangular  glass  microchannel,  which  is 
glued  to  a  circular  capillary  tube  and  a  170  pm  glass  coverslip  for  support.  Plastic 
tubing  connects  the  capillary  tube  to  the  syringe  pump. 


Figure  3.  Ensemble-averaged  velocity-vector  field  measured  in  a  30  pm  deep  x  300  pm  wide  x 
25  mm  channel.  The  spatial  resolution,  defined  by  the  interrogation  spot  size  of  the 
first  interrogation  window,  is  13.6  pm  x  4.4  pm  away  from  the  wall,  and  13.6  pm  x 
0.9  pm  near  the  wall.  A  5^%  overlap  between  interrogation  spots  yields  a  velocity 
vector  spacing  of  450  nm  in  the  wall-normal  direction  near  the  wall:  (a)  large-scale 
vector  field  showing  a  100  pm  x  85  pm  cross  section  of  the  channel,  (b)  near- wall 
view  of  the  lower  30  pm  of  vector  field. 
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Velocity  (mm/s) 

Figure  4.  Ensemble-averaged  velocity  profile  measured  in  a  nominally  30  pm  x  300  pm 
channel.  The  symbols  represent  ensemble-averaged  and  streamwise-averaged 
PIV  data.  The  solid  line  is  the  analytical  solution  for  Newtonian  flow  through 
a  rectangular  channel. 
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Figure  5.  Diagrams  depicting  the  different  ways  in  which  the  average  velocity  can  be  estimated:  (a) 
average  velocity  method,  (b)  average  image  method,  (c)  average  correlation  method. 


Figure  6.  Instantaneous  cross  correlation  functions  that  are  averaged  together  to  produce  an  average  correlatic 
function.  The  average  correlation  function  has  a  much  higher  signal-to-noise  ratio  than  tl 
instantaneous  correlation  functions. 
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Figure  7.  Velocity  vector  fields  showing  the  results  from  the  different  methods  of  calculation:  (a) 
Instantaneous  velocity  field,  (b)  time  average  of  twenty  instantaneous  velocity  fields,  (c) 
velocity  field  calculated  from  the  time  average  of  twenty  correlation  functions. 
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Figure  8.  Comparison  of  the  performance  of  the  three  averaging  techniques:  average  velocity  •, 
average  image  ■,  and  average  correlation  y. 
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Figure  9.  Comparison  of  a  forward  difference  adaptive  window  offset  analysis  (a)  versos  a 
central  difference  adaptive  window  offset  analysis  (b)  of  the  flow  around  a  single  red  blood  cell. 
This  figure  demonstrates  quantitatively  the  advantage  of  using  the  second  order  central 
difference  interrogation  (CDI)  scheme. 
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Figure  10.  Flowchart  demonstrating  iterative  procedure  necessary  adaptive  image  shifting  and 
second  order  accuracy. 
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Figure  11.  Schematic  views  illustrating  the  (a)  forward  difference  interrogation  algorithm  and  the 
(b)  central  difference  interrogation  algorithm  as  applied  to  the  cylindrical  Couette  flow  system. 
The  dashed  squares  represent  spatially  offset  interrogation  regions  while  the  dashed  circles 
represent  the  path  that  a  fluid  element  would  follow. 
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Camera 


Figure  12.  Experimental  setup.  The  camera-light  sheet  orientation  indicates  that 
velocity  components  in  a  radial-azimuthal  plane  will  be  measured. 
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Figure  13.  Convergence  characteristics  of  CDI  algorithm.  The  square  symbols  (■)  and  the  solid 
curve  show  the  RMS  difference  between  successive  iterations  normalized  by  the  RMS  of  the 
velocity  field.  The  circular  symbols  (•)  and  the  dashed  curve  show  the  fraction  of  the  velocity 
measurements  that  have  not  converged  as  a  function  of  iteration  number.  Both  data  sets  are 
modeled  quite  well  by  a  decaying  exponential  function. 
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Radial  position  (dimensionless) 


Figure  14.  Typical  distribution  of  velocity  measurements,  each  marked  by  a  plus  (+)  sign,  about 
the  predicted  bias  error  curve.  The  left  side  of  the  plot  corresponds  to  stationary  out^r  cylinder 
while  the  right  side  of  the  figure  corresponds  to  the  rotating  inner  cylinder.  The  scatter  of  the 
points  indicates  the  random  error  (±0.0048)  while  the  general  trend  indicates  the  velocity  bias 
error. 
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Inner  Cylinder  Rotation  (pixels) 


Figure  15.  Comparison  of  the  radial  velocity  errors  in  the  forward  difference  and  central 
difference  discretizations.  Square  symbols  (□,■)  represent  central  difference  interrogation  while 
circles  (O,#)  represent  forward  difference  interrogations.  Filled  symbols  (#,■)  represent 
interrogations  of  experimental  data  while  the  hollow  symbols  (0,D)  represent  the  Monte  Carlo 
results.  The  bottom  horizontal  axis  is  scaled  in  radians  of  rotation  of  the  inner  cylinder  while  the 
top  horizontal  axis  is  scaled  in  pixels  of  displacement  of  the  inner  cylinder.  The  left  vertical  axis 
is  scaled  as  the  fractional  velocity  error  while  the  right  vertical  axis  is  scaled  as  the  displacement 
error  in  pixels. 
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Displacement  Error  (pixels) 


Figure  16.  An  example  of  a  Silicon-micromachined  nozzle.  The  2-D  nozzle  contours  are 
machined  into  the  Silicon  using  Deep  Reactive  Ion  Etching  (DRIE).  The  top  and 
bottom  of  the  nozzle  are  formed  by  anodically  bonding  two  glass  wafers  to  the 
Silicon  wafer.  (Photo  courtesy  Robert  Bayt,  MIT) 
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Figure  17.  Schematic  of  a  micron-resolution  PIV  system.  The  system  incorporates  a  reflective 
DIC  technique  to  record  images  of  50  -  100  nm  dia.  particles.  An  optical  delay 
path  can  be  added  to  the  system  to  produce  two  pulses  of  laser  light  separated  by 
tens  of  nanoseconds. 
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(a)  (b) 


Figure  1 8.  Differential  Interference  Contrast  (DIC)  images  of  50  nm  dia.  Silicate  particles:  (a) 
recorded  using  a  NA  =  1.4  oil  immersion  lOOx  objective  lens,  (b)  recorded  using  a 
NA  =  0.6  air  immersion  40x  objective  lens 

Re  =  10,  3.4  Y  Resolution  96mm/sec-» 


Figure  19.  PIV  vector  field  of  water  flow  through  a  Silicon  micromachined  nozzle.  The  nozzle 
throat  is  70  pm  wide  throat  and  leads  into  a  constant-angle  expansion  section.  The 
spatial  resolution  of  the  PIV  measurements  are  10.9  x  3.4  pm  in  the  streamwise  and 
cross-flow  directions,  respectively. 
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